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Abstract 

This  report  presents  the  final  results  of  the  project  “Global  empirical  model  of  the  TBC  response  to 
geomagnetic  activity  and  forcing  from  below”  (FA8655- 12- 1-2057)  with  the  following  main 
objectives: 

a)  Development  of  global  empirical  background  TBC  model 

b)  Statistical  evaluation  of  the  global  background  TBC  model 

c)  Development  of  global  empirical  model  of  TBC  response  to  geomagnetic  activity 

d)  On-line  implementation  of  both  global  empirical  TBC  models:  background  and  TBC 
response  to  geomagnetic  activity  and  creating  of  their  mobile  versions. 

e)  Development  of  a  hybrid  global  TBC  model  -  an  attempt  to  mitigate  the  global  background 
TEC  model’s  error  by  using  regularly  arriving  new  CODE  TEC  data 

f)  Winter-time  assessment  of  the  global  TEC  dependence  on  the  stratospheric  temperature  and 
solar  radiation 


Introduction 

Description  of  the  problem 

The  high  sensitivity  of  the  ionosphere  to  the  external  forcing  related  to  changes  in  solar  EUV  and 
UV  radiation  and  geomagnetic  activity,  and  the  continuous  action  of  the  lower  atmospheric  forcing, 
as  the  wave  forcing  particularly  effective  during  low  solar  activity  conditions,  causes  significant 
ionospheric  variability  on  different  time  and  spatial  scales.  To  understand  and  forecast  such 
variability  is  the  main  tasks  of  space  weather  research.  Such  task  can  be  solved  by  building  of 
ionospheric  models  which  play  an  important  role  in  specifying  the  ionospheric  environment  as 
realistically  as  possible.  The  total  electron  content  (TEC)  has  received  a  great  deal  of  attention 
recently,  because  it  is  a  key  parameter  related  to  the  phase  delay  effects  on  the  signals  of  Global 
Navigation  Satellite  Systems  (GNSS).  The  ionospheric  effect  is  the  largest  error  source  in  GNSS 
positioning,  timing  and  navigation.  The  accuracy  of  the  GNSS  such  as  the  Global  Positioning 
System  (GPS),  the  Russian  GLONASS,  the  Chinese  BeiDou  and  the  European  Galileo,  is  heavily 
affected  by  the  ionosphere.  The  accurate  TEC  prediction,  particularly  during  periods  of  solar 
disturbances,  is  a  strong  requirement  for  the  reliable  performance  of  many  applications  including 
HF  communications,  satellite  positioning,  navigation  applications,  detection  and  tracking  of 
missiles  and  other  targets.  Most  of  the  ionospheric  error,  or  so  called  first-order  range  error,  has 
been  already  eliminated  by  differential  measurements  in  dual  frequency  systems  like  GPS,  1575.42 
MHz  at  LI  and  1227.60  MHz  at  L2.  However,  ionosphere  dual-frequency  algorithms  used  for 
positioning  applications  remove  only  first-order  range  error  but  do  not  take  into  account  its  higher-order 
terms.  Also,  the  ray  paths  and  TEC  are  assumed  to  be  the  same  for  both  frequencies  which  is  away 
from  the  reality  particularly  considering  the  horizontal  gradients  of  the  ionosphere  electron  density. 
Additionally,  there  are  still  numerous  single  frequency  applications  which  need  additional 
information  for  mitigating  the  ionospheric  propagation  error.  Such  GNSS  users  can  be  provided  by 
adequate  ionospheric  corrections  obtained  by  an  autonomous  ionospheric  TEC  model  (without  any 
ionospheric  measurements). 

Different  empirical  TEC  models,  based  on  existing  empirical  models  of  the  electron  density 
distribution  such  as  IRI  or  NeQuick  or  different  TEC  measurements,  have  been  already  built.  They 
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however  consider  only  the  external  forcing  of  the  ionosphere,  i.e.  variable  solar  BUY  and  UV 
radiation  and  geomagnetic  activity.  Often,  and  especially  recently  when  the  level  of  solar  activity  is 
very  low,  quite  large  day-to-day  changes  of  the  ionosphere  have  been  observed.  With  the  recent 
accumulation  of  satellite  measurements,  attention  is  now  being  directed  towards  investigating  the 
impact  of  the  processes  from  below  and  particularly  the  wave  forcing  from  the  lower  atmosphere. 
Various  ionospheric  observations  have  shown  the  development  of  longitudinal  wave-like  patterns 
(as  wave  number  three  (WN3),  or  wave  number  four  (WN4)  structures  that  revealed  quite  stable 
seasonal  behavior.  The  evidence  has  emerged  from  different  measurements  and  all  of  them 
unambiguously  display  manifestations  of  lower  atmospheric  dynamics  on  the  upper  atmosphere  and 
ionosphere.  Therefore,  it  is  time  a  new  type  of  global  empirical  TBC  model  to  be  built  where  both 
types  of  forcing,  from  above  and  from  below,  to  be  included. 


Summary  of  the  main  results 

In  order  to  build  a  new  type  of  global  empirical  TBC  model,  where  the  forcing  from  below  is 
included  as  well,  a  serious  knowledge  on  the  coupling  process  of  the  atmosphere-ionosphere  system 
is  requited.  The  generation  and  propagation  of  atmospheric  waves  is  a  dominant  aspect  of  the 
atmosphere  and  is  a  key  component  linking  different  regions.  Recent  studies  based  on  the 
observations  made  by  the  Sounding  of  the  Atmosphere  using  Broadband  Bmission  Radiometry 
(SABBR)  and  the  TIMBD  Doppler  Interferometer  (TIDI)  instruments  on  the  Thermosphere- 
lonosphere-Mesosphere-Bnergetics  and  Dynamics  (TIMBD)  satellite  have  provided  new  insight 
into  wave  fields  and  revealed  the  global  distribution  and  climatology  of  the  most  important  wave 
components  in  temperature  and  neutral  winds  respectively.  The  research  group  from  the  NIGGG- 
BAS  has  serious  contribution  for  clarifying  the  global  distribution  and  temporal  variability  of  the 
main  atmospheric  wave  fields,  as  tides  and  planetary  waves,  which  have  impact  on  the  ionosphere. 
Pancheva  et  al.  (2009a,  2009b)  created  an  advanced  method  for  analysis  of  satellite  measurements 
on  the  basis  of  which  the  authors  have  been  able  to  study  in  detail  the  wave  forcing  of  the 
ionosphere  from  lower  stratosphere  to  lower  thermosphere,  i.e.  up  to  the  dynamo  region  where  the 
waves  directly  or  indirectly  (by  electrodynamics)  can  have  impact  on  the  ionosphere.  By  using  this 
advanced  approach  the  climatological  features  of  the  following  atmospheric  waves  have  been 
investigated:  migrating  diurnal  tide  (Mukhtarov  et  al.,  2009),  migrating  semidiurnal  tide  (Pancheva 
et  al.,  2009c),  nonmigrating  tides  (Pancheva  et  al.,  2009d;  2010a)  and  different  stationary  and 
zonally  propagating  planetary  waves  (Pancheva  et  al.,  2007,  2009b,  2010b;  Mukhtarov  et  al., 
2010a;  Lu  et  al.,  2012).  Very  recently  the  climatology  of  the  migrating  terdiurnal  tide  has  been 
reported  also  (Pancheva  et  al.,  2013).  The  detailed  picture  of  the  spatial  structure  and  temporal 
variability  of  the  atmospheric  tides  and  planetary  waves  obtained  from  satellite  measurements  and 
summarized  in  the  Springer  book  chapter  by  Pancheva  and  Mukhtarov  (2011a)  has  already  served 
as  a  benchmark  and  guide  for  numerical  modeling  studies  aimed  at  better  understanding  the 
coupling  processes  by  tidal  and  planetary  wave  patterns. 

The  next  step  in  the  atmosphere-ionosphere  coupling  studies  is  clarifying  the  ionospheric  response 
to  different  tides  and  planetary  waves.  For  this  purpose  the  ionospheric  measurements  of  the 
Constellation  Observing  System  for  Meteorology,  Ionosphere,  and  Climate  (COSMIC),  as  the 
COSMIC  electron  density  profiles  in  the  altitude  range  of  100-800  km,  have  been  used.  The  crucial 
point  in  studying  the  ionospheric  response  to  wave  forcing  from  below  is  the  application  of  one  and 
the  same  data  analysis  method  to  both  data  sets,  atmospheric  (SABBR/TEMBD  satellite  data)  and 
ionospheric  (COSMIC  electron  density)  data.  By  using  two  data  sets  Pancheva  and  Mukhtarov 
(2010)  for  the  first  time  provided  evidence  showing  that  the  ionospheric  WN4  and  WN3  are 
generated  mainly  by  eastward  propagating  nonmigrating  diurnal  tides  with  wavenumbers  3  and  2 
respectively.  The  spatial  structures  of  the  ionospheric  response  to  some  migrating  and  nonmigrating 
tides  have  been  presented  in  Mukhtarov  and  Pancheva  (2011)  and  Pancheva  and  Mukhtarov 
(2012a),  while  the  ionospheric  response  to  planetary  waves  is  considered  in  Mukhtarov  et  al. 
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(2010b)  and  Pancheva  and  Mukhtarov  (2012b).  The  detailed  picture  of  the  spatial  structure  and 
temporal  variability  of  the  ionospheric  response  to  waves  forced  from  below  has  been  summarized 
in  the  Springer  book  chapter  by  Pancheva  and  Mukhtarov  (2012c).  The  obtained  main  features  of 
the  ionospheric  response  to  different  tides  forced  from  below  have  been  supported  by  the  first 
detailed  comparison  between  simulated,  by  the  whole  atmospheric  model  called  GAIA,  and 
observed  COSMIC  global  electron  density  responses  (Pancheva  et  al.,  2012).  Recently  a  special 
attention  has  been  paid  to  the  vertical  coupling  during  anomalous  stratospheric  events  like  sudden 
stratospheric  warmings  (SSW).  Pancheva  and  Mukhtarov  (2011b)  by  using  atmospheric  and 
ionospheric  satellite  measurements  for  the  first  time  presented  experimental  evidence  that  the 
ionosphere  regularly  responds  to  almost  all  SSW  stratospheric  temperature  pulses  at  high 
stratospheric  latitudes.  Later  by  combining  the  observations  with  simulations  done  by  the  GAIA 
model  further  insight  on  the  ionospheric  response  to  the  SSW  events  has  been  obtained  (Jin  et  al., 
2013). 

Besides  knowledge  on  atmosphere-ionosphere  coupling  the  NIGGG-BAS  research  group  has 
serious  experience  in  empirical  modeling  also.  It  participated  in  almost  all  COST  actions  devoted  to 
the  applied  ionospheric  research  and  specialized  in  empirical  modeling  of  ionosphere  and 
development  of  methods  for  nowcasting  and  forecasting  ionospheric  parameters  of  interest  to  radio 
communications  and  GNSS  navigation.  The  single-station  model  approach  for  long-term  prediction 
was  developed  quite  long  time  ago  by  Pancheva  and  Mukhtarov  (1996,  1998).  The  idea  for 
describing  the  solar  activity  by  two  parameters,  i.e.  the  level  of  solar  activity  and  its  tendency,  was 
introduced  for  the  first  time  by  Pancheva  and  Mukhtarov  (1996)  in  modeling  the  monthly  median 
critical  frequency  of  the  ionospheric  F-region,/oF2,  above  Sofia  and  is  used  in  the  newly  developed 
background  TBC  model  (Section  1).  A  method  for  “weighted  extrapolation”  of  past  measurements 
of  the  foF2  has  been  developed  by  using  the  linear  regression  method,  which  coefficients  were 
obtained  from  its  autocorrelation  function  (Muhtarov  and  Kutiev,  1999).  Based  on  this  method,  a 
series  of  single-station  models  for  short-term  forecasting  (up  to  3  days  ahead)  of  foF2  have  been 
developed.  The  basic  model  (Kutiev  et  al,  1999,  Muhtarov  et  al,  2001)  has  been  implemented  in 
STIF  software  at  Rutherford  Appleton  Laboratory,  UK.  Muhtarov  et  al,  (2001)  introduced  a  new 
term  in  the  autoregressive  formula  for  prediction  of  foF2,  representing  the  changes  of  geomagnetic 
activity.  This,  so  called  geomagnetically-correlated  autoregression  method  was  implemented  in  the 
forecasting  software  of  DIAS  system  (http://www.iono.noa.gr/DIAS).  Important  improvement  of 
forecasting  methodology  was  made  by  separation  of  the  seasonal  variations  of  foF2,  represented  by 
monthly  medians,  from  the  deviations  presumably  invoked  by  geomagnetic  activity;  this  approach 
is  applied  in  the  created  new  global  TEC  model  response  to  geomagnetic  activity  (Section  3).  The 
relative  deviations  of  foF2  (denoted  as  rfoF2)  from  its  median  values  were  represented  by  analytical 
functions  of  the  geomagnetic  Kp-index.  Another  models,  as  the  midlatitude  model  of  foF2  (Kutiev 
and  Muhtarov,  2001)  and  the  global  model  of  monthly  average  deviations  as  function  of  Kp  (Kutiev 
and  Muhtarov,  2003)  have  been  developed  in  the  framework  of  COST  actions  and  reflected  the 
European  level  of  applied  ionospheric  research,  which  in  many  cases  was  leading  in  the  global 
scale. 

All  above  mentioned  studies  have  been  related  mainly  to  the  ionospheric  parameters  defined  by 
ionosonde  stations.  Recently  Andonov  et  al.  (2011)  have  presented  an  empirical  TEC  model 
response  to  the  geomagnetic  activity  for  the  American  sector  but  particularly  valid  for  low  solar 
activity.  It  was  based  on  the  two-dimensional  (2D)  cross-correlation  analysis  which  revealed  both 
positive  and  negative  phases  of  the  response.  The  both  phases  of  the  ionospheric  response  have 
different  duration  and  time  delay  with  respect  to  the  geomagnetic  activity,  season  and  geographical 
latitude.  The  same  approach  is  applied  and  developed  further  in  the  created  new  global  TEC  model 
response  to  geomagnetic  activity  (Section  3). 

Based  on  the  above  described  knowledge  and  experience  two  global  empirical  TEC  models  have 
been  built  during  the  first  year  of  the  project.  The  first  model  is  a  background  TEC  model 
(Mukhtarov  et  al.,  2013a)  for  describing  the  mean  behavior  of  the  ionosphere  under  both  its 
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primary  external  driver,  i.e.  the  direct  photo-ionization  by  incident  solar  radiation,  and  regular  wave 
particularly  tidal  forcing  from  the  lower  atmosphere.  Moreover,  we  aim  to  make  this  model 
applicable  under  quiet  geomagnetic  conditions  for  long-term  prediction  of  the  average  TBC 
variability.  This  model  is  statistically  evaluated  and  even  an  error  model  is  suggested  (Mukhtarov  et 
al.,  2013b).  The  second  model  is  a  global  empirical  model  of  the  TBC  response  to  geomagnetic 
activity  (Mukhtarov  et  al.,  2013c)  described  by  the  Kp-index  which  is  aimed  at  short-term 
prediction  (a  day  ahead)  of  the  TBC  variability. 


1.  Development  of  empirical  background  TEC  model 


1.1  TEC  Data  Set 


It  is  known  that  empirical  models  typically  represent  the  gross  features  in  the  ionosphere  quite  well, 
but  are  limited  to  the  way  the  model  was  built,  the  data  that  was  used  to  construct  it,  and  the 
conditions  that  were  occurring  while  the  data  was  taken.  Recently,  the  GPS  measurements  obtained 
from  the  global  and  regional  networks  of  International  GNSS  Service  (IGS)  ground  receivers  have 
become  a  major  source  of  TBC  data  over  large  geographic  areas.  This  system  offers  low  cost 
information  characterized  by  its  accuracy,  high  temporal  and  spatial  resolution,  and  availability. 
The  present  background  TBC  model  is  constructed  on  the  base  of  vertical  TBC  maps  generated  by 
the  CODB  at  Astronomical  and  Physical  Institutes  of  the  University  of  Bern 
(http://cmslive3.unibe.ch/unibe/philnat/aiub/content/el5/e59/el26/e440/e447/index_eng.html).  We 
particularly  note  that  TBC  everywhere  means  vertical  TEC.  The  data  for  full  13  years,  1  January 
1999  -  31  December  2011,  provided  from  the  CODE  FTP  directory:  ftp://ftp.unibe.ch/aiub/CODE/ 
are  used.  The  two-hourly  sets  are  derived  from  GPS  data  of  the  global  IGS  network  of  about  200 
stations.  The  GIM/CODE  is  regarded  as  one  of  the  most  precise  TEC  maps  generated  from  GNSS 
observations.  The  used  global  IGS  TEC  data  have  a  time  resolution  of  2  h  and  a  grid  spacing  of  5° 
X  2.5°  in  longitude  and  latitude,  respectively  with  errors  of  several  TEC  Units  (TECU,  1  TECU  = 
10^^  el/m^). 


The  original  global  TEC  data  were  arrayed  in  terms  of  the  coordinate  system  of  geographical 
latitude  (from  -87.5°  to  87.5°  at  each  2.5°)  and  longitude  (from  -180°  to  180°  at  each  5°).  It  is  known 
however  that  the  neutral  wind  and  electric  field  effects  on  the  ionosphere  are  dependent  on  the 
geomagnetic  field  configuration  as  the  electrons  are  constrained  to  the  magnetic  field  lines.  That  is 
why  the  distribution  of  the  ionospheric  parameters,  including  TEC  as  well,  is  usually  presented  in 
geomagnetic  latitude  instead  of  geographic  one.  We  use  the  modified  dip  latitude  (modip), 
introduced  by  Rawer  (1963).  The  modified  dip  (modip)  latitude  which  is  adapted  to  the  real 


magnetic  field,  e.g.,  to  the  magnetic  inclination  (dip),  is  defined  by:  tan//  = 


,  where  is 


the  modip  latitude,  1  is  the  true  magnetic  dip  (usually  at  a  height  of  350  km),  and  (P  is  the 
geographic  latitude.  Modip  equator  is  the  locus  of  points  where  the  magnetic  dip  (or  inclination)  is 
0.  In  the  equatorial  zone,  the  lines  of  constant  modip  are  practically  identical  to  those  of  the 
magnetic  inclination  but  as  latitude  increases  they  deviate  and  come  nearer  to  those  of  constant 
geographical  latitude.  The  poles  are  identical  to  the  geographic  ones.  For  the  purpose  of  this  study 
the  global  TEC  data  were  re-arrayed  in  terms  of  the  coordinate  system  of  modip  latitude,  from  -80° 
to  80°  at  each  5°,  and  geographic  longitude,  from  -180°  to  180°  at  each  15°.  The  TEC  data  falling 
into  the  area  5°  (modip  latitude)  x  15°  (longitude)  were  averaged. 


Usually  the  background  ionospheric  models  are  formulated  in  terms  of  monthly  median  parameters 
because  such  parameters  are  not  affected  by  large  but  short-time  lasting  disturbances  generated  by 
strong  geomagnetic  storms.  In  the  present  study  we  use  sliding  medians  defined  by  a  31 -day 
moving  window  and  the  median  value  is  assigned  to  the  central  day  of  the  window,  i.e.  to  the  16* 
day  of  the  window.  The  sliding  medians  are  calculated  independently  for  each  point  of  the  grid  (as 
it  is  done  with  single  station  data).  In  this  way  the  daily  TEC  time  series  are  obtained  at  each  modip 
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latitude,  longitude  and  UT.  It  is  worth  noting  that  the  perturbations  from  geomagnetic  origin  or 
related  to  solar  rotation  period  are  filtered  from  these  time  series. 

1.2  Basic  Approach  of  the  Model  Construction 

The  basic  idea  of  each  global  empirical  background  TBC  model,  used  for  long-term  prediction,  is  to 
define  a  set  of  empirical  functions  which  describes  the  most  probable  TBC  values  at  given  solar 
activity,  day  of  the  year  (DOY),  geographic  location  and  UT.  In  the  present  study  we  accepted:  (i) 
longitude  and  UT  as  independent  variable  quantities;  the  conversion  into  FT  is  a  simple  procedure, 
and  (ii)  at  each  modip  latitude  a  separate  model  is  constructed;  the  values  of  the  model  TBC  which 
do  not  belong  to  the  5°  modip  grid  are  obtained  by  an  interpolation  procedure  that  will  be  described 
later.  The  latter  is  done  because  if  a  latitudinal  approximation  is  used  first  the  number  of  model 
constants  will  increase  and  second  an  additional  error  will  be  introduced  in  the  model. 

According  to  the  above  mentioned  approach  the  TBC  values  at  each  modip  latitude  circle  can  be 
presented  as  a  function  of: 

TEC  (solar  activity,  DOY,  longitude,  UT)  (1) 

Ideally  the  solar  activity  should  be  described  by  an  index  that  tracks  the  solar  cycle  changes  in  the 
BUY  wavelength  range,  since  this  part  of  the  solar  spectrum  affects  the  ionosphere.  However,  such 
indices  cannot  be  observed  at  the  ground  and  are  only  available  for  relatively  short  time  periods 
covered  by  satellite  UV  instruments.  Thus,  most  ionospheric  modelers  use  the  sunspot  number 
(number  of  dark  spots  on  the  solar  disc)  and  the  solar  radio  flux  at  10.7  cm  wavelength  (F10.7)  as 
solar  indices,  since  both  can  be  observed  from  the  ground,  long  data  records  exist  and  they  can  be 
predicted.  These  indices  together  with  their  6-month  predictions  are  regularly  published  by  NOAA 
Space  Weather  Prediction  Center  (http://www.swpc.noaa.gov/index.html).  In  the  present  study 
F10.7  is  used  as  a  proxy  for  the  solar  activity.  It  is  known  however  that  the  ionosphere  behaves 
differently  at  the  rising  and  declining  phase  of  the  solar  cycle  at  one  and  the  same  F10.7.  To  include 
this  ionospheric  feature  in  the  model  an  additional  parameter  Kf  is  used  which  describes  the  linear 
rate  of  change  of  FI 0.7.  As  has  been  already  mentioned  this  idea  was  introduced  for  the  first  time 
by  Pancheva  and  Mukhtarov  (1996).  Figure  1  shows  the  temporal  variability  of  the  used  two  solar 
parameters  FI  0.7  in  solar  flux  units  (lO'^^W  m'^  Hz'^)  and  Kp  for  the  considered  13  years  (1999- 
2011). 


Figure  1  Monthly  mean  solar  radio  flux  F10.7  (red  line)  and  its  linear  rate  of  change  Kp  (blue  line)  obtained  by  a 
sliding  11 -month  time  derivatives  of  FI 0.7. 


The  seasonal  variability  of  the  vertical  TBC  is  characterized  mainly  by  annual  and  semiannual 
components  whose  amplitudes  and  phases  depend  on  solar  activity,  geographic/geomagnetic 
location  and  UT.  The  seasonal  components  with  periods  shorter  than  6  months  have  also  some 
contribution  but  they  are  weaker  than  annual  and  semiannual  components.  In  general,  the 
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semiannual  behavior  is  characterized  by  a  larger  peak  in  March-April  than  that  in  September- 
October;  the  difference  is  particularly  strong  (up  to  40%)  during  high  solar  activity. 

In  all  existing  so  far  empirical  models  the  diurnal  TBC  variability  is  described  only  by  the  migrating 
“tidal”  components,  i.e.  components  with  periods  24,  12,  8  hours,  etc.,  which  propagate  with  the 
apparent  motion  of  the  Sun  to  a  ground-based  observer.  In  this  way  it  is  assumed  that  the  diurnal 
variability  of  the  ionosphere  is  forced  mainly  by  the  diurnal  cycle  of  the  photo-ionization.  However, 
only  within  the  past  5-6  years  has  the  realization  emerged  that  “troposphere  weather”  contributes 
significantly  to  the  “space  weather”  of  the  thermosphere,  especially  during  solar  minimum 
conditions.  Much  of  the  attendant  variability  is  attributable  to  upward-propagating  solar  tides 
excited  by  latent  heating  due  to  deep  tropical  convection,  and  solar  radiation  absorption  primarily 
by  water  vapor  and  ozone  in  the  troposphere  and  stratosphere,  respectively.  Recent  studies  based  on 
the  modern  satellite-board  data  (electron  densities  from  COSMIC  and  temperatures  from 
SABBR/TEMBD)  revealed  the  importance  of  the  ionospheric  forcing  from  below  (see  previous 
section).  There  is  another  reason  also  for  including  the  nonmigrating  tidal  variability  in  the  TBC 
model.  It  is  related  to  the  offset  between  the  geographic  and  modip  latitudes.  The  dynamics  of  the 
thermosphere  (defined  mainly  by  prevailing  winds  and  atmospheric  tides)  is  defined  in  a  geographic 
frame  but  its  effect  on  the  ionosphere  depends  on  the  geomagnetic  field  configuration.  The  photo¬ 
ionization  depends  also  on  geographic  frame  (solar  zenith  angle)  however  as  the  electrons  are 
constrained  to  the  magnetic  field  lines  its  effect  is  geomagnetic  field  dependent.  The  offset  between 
geographic  and  geomagnetic  frames  generate  additional  tidal  ionospheric  responses  which  can 
originally  not  be  present  in  the  forcing  from  below.  These  additional  ionospheric  tides  are  much 
weaker  than  the  real  ones  (on  the  average  ~3-10%)  but  they  have  some  contribution  to  the  diurnal 
variability  of  the  ionosphere  particularly  in  shaping  some  quasi-stationary  structures. 

As  the  time  scales  of  the  solar  cycle,  seasonal  and  diurnal  influences  on  the  TBC  variability  are  very 
different  (they  differ  at  least  an  order  of  magnitude)  then  the  shorter-period  TBC  variabilities  are 
usually  modulated  by  the  longer  ones.  In  this  case  the  TBC  spatial-temporal  variability  can  be 
represented  as  a  multiplication  of  three  separable  functions,  i.e.  (1)  can  be  expressed  in  the 
following  way: 

TEC{F  ,Kp,DOY  ,lon,UT)  =  ^,{P  ^KF)^2i^OY)^^{lon,UT)  (2) 


The  above  right  hand  side  unknown  functions  0k  (k=l,2,3)  can  be  represented  by  their  series 
expansions;  0]  can  be  expanded  in  Taylor  series,  while  02  and  03,  which  are  periodic  functions 
with  periods  respectively  a  year  and  a  solar  day,  can  be  expanded  in  Fourier  series.  Therefore,  the 
background  TBC  model  can  be  described  by  the  following  functions: 
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The  expression  in  the  first  right  hand  bracket,  i.e.  the  Taylor  series  expansion  up  to  degree  of  2, 
represents  the  solar  activity  term  which  modulates  the  seasonal  and  diurnal  behavior  of  the 
ionosphere.  The  seasonal  term  (expression  in  the  second  right  hand  bracket)  includes  4 
subharmonics  of  the  year,  i.e.  annual,  semiannual,  4-  and  3-month  components;  it  modulates  the 
diurnal  behavior  of  the  ionosphere.  The  diurnal  variability  of  the  TBC  model  (expression  in  the 
third  right  hand  bracket)  is  composed  by  three  parts:  mean  TBC  (yo),  a  part  describing  migrating 
and  nonmigrating  tides  and  a  part  representing  the  effect  of  the  stationary  planetary  waves  (SPWs). 
The  contribution  of  the  migrating  and  nonmigrating  tides  is  presented  by  2D  (longitude-time)  sine 
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functions  with  zonal  wavenumbers  up  to  4  and  4  subharmonics  of  the  solar  day  (24-,  12-,  8-  and  6- 
hour  components).  The  last  part,  describing  the  contribution  of  the  SPWs,  includes  waves  with 
zonal  wavenumbers  up  to  4.  The  presence  of  these  wave  structures  in  the  ionosphere  can  be  related 
to  a  few  reasons:  (i)  offset  between  geographic  and  modip  frames;  (ii)  can  be  generated  by  coupling 
processes  between  migrating  and  nonmigrating  tides  with  one  and  the  same  periods,  and  (iii)  some 
effect  of  the  SPW 1  temperature  field  in  the  lower  thermosphere  on  the  ionosphere,  particularly  at 
middle-high  and  high  latitude  ionosphere;  Mukhtarov  et  al.  (2010a)  found  strong  evidence 
indicating  that  the  auroral  heating  is  a  main  origin  of  the  lower  thermospheric  SPW  1  structure. 


The  background  TBC  model  described  by  (3)  contains  4374  constants  and  they  are  determined  by 
least  squares  fitting  techniques.  The  numbers  of  the  included  components  in  the  above  described 
Taylor  and  Fourier  expansion  series  are  defined  experimentally.  We  accepted  only  the  above 
mentioned  solar,  seasonal  and  diurnal  components  because  the  addition  of  more  components  does 
not  improve  the  model  error. 


Month  Number  (start  January  1999) 


Month  Number  (start  January  1999) 


Month  Number  (start  January  1999) 


Figure  2  Latitude-time  cross-sections  of  the  following  TEC  components:  (i)  zonal  and  time  mean  TEC  (upper  most 
plot);  (ii)  left  column  of  plots  shows  the  amplitudes  of  migrating  diurnal  (DWl,  upper  plot),  semidiurnal  (SW2,  middle 
plot)  and  terdiurnal  (TW3,  bottom  plot)  tides,  and  (iii)  right  column  of  plots  shows  the  amplitudes  of  SPWl  (upper  plot), 
nonmigrating  zonally  symmetric  diurnal  tide  (DO,  middle  plot)  and  nonmigrating  eastward  propagating  tide  with  zonal 
wavenumber  3  (DE3,  bottom  plot).  The  considered  years  from  1999  to  2011  are  separated  by  thin  white  lines. 
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Figure  2  presents  examples  of  solar  cycle  and  seasonal  modulations  of  some  diurnal  components 
included  in  the  TBC  model.  For  this  purpose  the  TBC  data  are  decomposed  to  mean  TBC,  migrating 
and  nonmigrating  tides,  and  SPWs  (i.e.  expressions  in  the  third  right  hand  bracket  of  (3))  by  using  a 
31 -day  window.  Then  the  31 -day  window  is  moved  forward  through  the  time  series  with  steps  of  1 
day  in  order  to  obtain  the  daily  values  of  the  wave  characteristics  for  period  of  time  1  January  1999- 
31  December  2011.  The  monthly  mean  wave  characteristics,  shown  in  Figure  2,  are  obtained  by 
vector  averaging  for  each  calendar  month.  The  uppermost  plot  shows  the  latitude-time  cross-section 
of  the  zonal  and  time  mean  of  the  TBC  (first  term,  yo,  in  the  third  bracket);  the  considered  years 
from  1999  to  2011  are  separated  by  thin  white  lines.  This  diurnal  TBC  component  follows  strictly 
the  solar  activity,  even  the  two  maxima,  a  main  maximum  in  2002  (~80  TBCU)  and  a  secondary 
one  in  2000  (~75  TBCU),  can  be  clearly  distinguished.  The  semiannual  variability  is  a  dominant 
component  of  the  seasonal  behavior;  on  the  average  the  vernal  equinox  maxima  are  stronger  than 
the  autumnal  ones.  The  left  column  of  plots  present  the  latitude-time  cross-sections  of  the  first  three 
migrating  tides,  24-h  (DWl,  upper  plot),  12-h  (SW2,  middle  plot)  and  8-h  (TW3,  bottom  plot).  The 
DW 1  component  is  the  strongest  tidal  component  (maximum  amplitude  reaches  ~50  TBCU)  and  is 
shaped  mainly  by  the  diurnal  variability  of  the  photo-ionization.  The  SW2  component  is 
significantly  weaker  than  the  DWl  one  (maximum  amplitude  reaches  ~11  TBCU)  and  as  it  has  been 
already  mentioned  is  formed  mainly  by  the  SW2  tide  forced  from  the  lower  atmosphere  (Pancheva 
and  Mukhtarov,  2012a).  A  clear  evidence  for  the  dominant  role  of  the  lower  atmospheric  SW2  tide 
on  the  TBC  SW2  variability  is  the  existence  of  the  local  winter  maxima  at  modip  latitude  around 
±60“.  The  solar  cycle  and  seasonal  (mainly  semiannual)  modulations  can  be  seen  well  for  these 
diurnal  components  as  well.  While  at  high  solar  activity  the  separation  at  both  sides  of  the  equator 
can  be  well  seen  for  the  mean  (yo)  and  DWl  components  for  the  SW2  component  such  separation  is 
evident  better  at  low  solar  activity;  at  high  solar  activity  the  SW2  amplifies  predominantly  over  the 
equator.  The  amplifications  of  the  mean  TEC,  tidal  DWl  and  SW2  amplitudes  around  ±(20-30“) 
modip  is  related  to  the  equatorial  ionization  anomaly  (El A)  observed  mainly  during  the  daytime. 
Figure  2  shows  that  the  solar  activity  affects  not  only  the  amplitude  of  the  TEC  equatorial  anomaly 
but  also  the  location  of  the  crests;  at  high  solar  activity  they  are  located  close  to  ±30°  while  at  low 
solar  activity  the  crests  move  close  to  the  equator,  around  ±20°.  The  solar  cycle  and  seasonal 
modulations  are  seen  also  on  the  third  migrating  component,  TW3.  The  latitude  structure  of  this 
diurnal  component  shows  a  main  amplification  over  the  equator  and  secondary  ones  around  ±50° 
and  ±65°;  the  latter  are  particularly  well  seen  during  high  solar  activity. 

The  right  column  of  plots  present  the  latitude-time  cross-sections  of  the  amplitudes  of  the  SPWl 
(upper  plot),  nonmigrating  zonally  symmetric  diurnal  (DO,  middle  plot)  and  nonmigrating  eastward 
propagating  diurnal  tide  with  zonal  wavenumber  3  (DE3,  bottom  plot).  All  these  diurnal 
components  demonstrate  regular  solar  cycle  and  seasonal  variability.  Their  amplitudes  are  weaker 
than  the  DWl  tide  but  are  comparable  with  those  of  the  SW2  and  TW3  tides.  The  DO  is  the 
strongest  nonmigrating  component  with  maximum  amplitude  of  ~11  TBCU  (the  same  as  that  of 
SW2).  It  amplifies  mainly  in  the  Southern  Hemisphere  (SH)  at  high  latitudes;  similar  distributions 
have  also  the  other  zonally  symmetric  tidal  components  but  their  amplitudes  are  weaker  than  that  of 
DO  (not  shown  here).  All  zonally  symmetric  tidal  components  show  amplifications  like  stripes 
between  -40°  and  -70°  modip  latitude  that  can  be  distinguished  even  at  low  solar  activity.  Similar 
amplifications  are  evident  for  SPWs  components,  but  they  are  present  at  both  hemispheres.  Later  it 
will  be  shown  that  just  these  zonally  symmetric  nonmigrating  and  SPW  components  have 
predominant  contribution  to  the  so  called  Weddell  Sea  Anomaly  (WSA).  This  anomaly  appears  as 
an  evening  enhancement  in  electron  density,  i.e.  larger  nighttime  electron  density  than  during  the 
day,  in  the  region  near  the  Weddell  Sea,  Antarctica  peninsula.  The  WSA  occurs  mostly  in  southern 
summer  and  can  extend  from  South  America  and  Antarctica  to  the  central  Pacific.  The  TEC  DE3 
component  (bottom  plot),  similarly  to  the  other  diurnal  components,  is  strongly  modulated  by  the 
solar  cycle  and  shows  clear  semiannual  variability.  In  this  case  however  the  autumnal  maxima  are 
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stronger  than  the  vernal  ones.  This  is  due  to  fact  that  the  neutral  DB3  tide  forced  from  below  is  the 


main  driver  of  the  ionospheric  DBS  variability  (Pancheva  and  Mukhtarov,  2010);  the  authors 
presented  unambiguous  evidence  that  the  ionospheric  DBS  variability  is  forced  from  below.  We  pay 
special  attention  to  the  TBC  DBS  nonmigrating  diurnal  component  as  it  has  the  main  contribution  to 
the  so  called  WN4  longitude  structure  observed  in  many  ionospheric  parameters  particularly  during 
low  solar  activity.  Later  we  will  demonstrate  the  ability  of  this  background  model  to  reproduce  the 
WN4  and  WNS  ionospheric  structures. 

1.3  Model  Results 

The  basic  aim  of  each  global  TBC  model  used  for  long-term  prediction  is  to  construct  the  global 
distribution  of  the  TBC,  i.e.  to  obtain  global  TEC  maps,  at  given  solar  activity,  day  of  the  year  and 
UT.  The  TEC  maps  are  constructed  by  interpolation  of  the  TEC  values  from  the  used  grid  with  a  5° 
step  in  modip  latitude  and  15°  in  longitude.  The  interpolation  between  obtained  from  the  model 
TEC  values  is  done  by  using  Inverse  Distance  Method  {Shepard,  D.,  1968).  Then  the  modip  frame 
is  converted  to  geographical  one.  The  TEC  values  assigned  to  both  poles  are  found  by  interpolation 
between  the  known  from  the  model  points  which  have  the  highest  northern  and  southern  latitudes. 
The  model  maps  are  arrayed  in  terms  of  the  coordinate  system  of  geographical  latitude  from  -90°  to 
90°  at  each  5°  and  longitude  from  -180°  to  180°  at  each  5°. 

First  we  will  show  how  the  background  TEC  model  describes  the  WSA.  The  zone  of  anomalous 
diurnal  variations  in  foF2,  which  is  characterized  by  an  excess  of  nighttime  foF2  values  over 
daytime  ones,  occupies  the  longitudes  of  0°-180°W  and  the  latitudes  of  40°-80°S  as  the  effect  is 
maximal  (up  to  ~5  MHz  for  the  critical  frequency  of  the  F-region, /0F2)  at  longitudes  of  40°-105°W 
and  latitudes  of  60°-70°S.  Figure  S  presents  the  global  map  in  geographical  coordinate  system 
calculated  from  the  empirical  background  TEC  model  (upper  plot)  and  compared  with  the  CODE 
TEC  map  (bottom  plot)  at  08UT  for  12  December  2012.  The  stripe  TEC  amplification  in  the 
Western  Hemisphere  at  latitudes  of  ~50-80°S,  i.e.  the  WSA,  can  be  clearly  distinguished  at  both 
model  and  CODE  TEC  maps;  the  maximal  effect  at  both  maps  is  near  70°S  and  longitudes  of  ~0- 
100°W.  The  presence  of  the  WSA  is  a  reason  for  appearing  of  an  additional  to  the  equatorial 
anomaly  TEC  amplification  around  S0-40°S  and  at  the  most  Western  Hemisphere;  this  feature  is 
also  well  reproduced  by  the  model.  The  model  TEC  map  describes  well  also  the  hemispheric 
asymmetry  of  the  equatorial  anomaly  revealing  that  the  summer  crest  is  stronger  than  the  winter  one. 
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Figure  3  Global  map  in  geographical  frame  calculated  from  the  empirical  background  TEC  model  (upper  plot) 


compared  with  the  CODE  TEC  map  (bottom  plot)  at  08  UT  for  12  December  2012. 

To  demonstrate  how  the  model  reproduces  some  longitude  wave-like  structures  we  re-arrayed  the 


model  and  CODE  TEC  data  sets  in  local  time  (LT).  It  has  been  already  mentioned  that  usually  the 


ionospheric  wave-like  longitude  structures  are  observed  during  low  solar  activity  as  the  WN4  is 
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seen  in  August-October  while  the  WN3  in  December- January.  Figure  4  shows  the  comparison 
between  the  global  TEC  model  maps  in  modip  latitude  (upper  row  of  plots)  and  CODE  TEC  maps 
(bottom  row  of  plots)  which  represent  the  ionospheric  WN4  (left  column  of  plots)  and  WN3  (right 
column  of  plots)  structures.  The  example  for  the  WN4  structure  is  for  October  2008  at  23ET,  i.e. 
night-time,  and  that  is  why  there  is  no  the  splitting  of  the  irregularities  at  both  sides  of  the  equator. 
Four  peaks  around  longitudes:  -150°,  -60°,  30°  and  120°  can  be  clearly  distinguished  at  both  model 
and  CODE  TEC  maps.  There  is  not  only  qualitative  but  also  quantitative  similarity  between  the 
model  and  observations.  The  example  for  the  WN3  structure  is  for  January  2008  at  14LT,  i.e. 
daytime,  when  the  equatorial  anomaly  is  present.  The  splitting  of  the  irregularities  at  both  sides  of 
the  equator  is  seen  at  both  model  and  CODE  TEC  maps.  However,  the  WN3  structure  is  well 
developed  and  significantly  stronger  in  SH.  This  hemispheric  asymmetry  is  due  to  the  asymmetry  of 
the  ionospheric  DE2  variability,  reported  by  Mukhtarov  and  Pancheva  (2011),  which  is  the  main 
contributor  of  the  WN3  structure.  The  three  peaks  particularly  in  the  SH  are  located  at  longitudes  of 
-150°,  0°  and  120°.  They  are  not  exactly  equidistant  most  probably  because  the  contribution  of  other 
nonmigrating  tides,  as  DW4  and  SEl,  and  SPW3  (Pancheva  and  Mukhtarov,  2012a,  2012c).  There 
are  some  signatures  for  the  first  and  third  peaks  at  NH  evident  at  both  model  and  CODE  TEC  maps. 
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Figure  4  Comparison  between  the  global  TEC  model  maps  in  modip  latitude  (upper  row  of  plots)  and  CODE  TEC 
maps  (bottom  row  of  plots)  which  shows  the  ionospheric  WN4  (left  column  of  plots,  at  23  LT  for  15  October  2008)  and 
WN3  (right  column  of  plots,  at  14  LT  for  15  January  2008)  structures. 


The  ability  of  the  background  TEC  model  to  reproduce  the  temporal-spatial  feature  of  the  input 
CODE  TEC  data  will  be  displayed  by  a  comparison  between  the  model  and  CODE  TEC  maps  for 
different  solar  activity,  seasons  and  UT.  While  Figure  5a  shows  the  global  maps  in  geographical 
coordinate  system  calculated  from  the  empirical  background  TEC  model  (left  column  of  plots) 
which  are  compared  with  the  CODE  TEC  maps  (right  column  of  plots)  at  12UT  for  15  January 
2001  (upper  row  of  plots)  and  15  March  2001  (bottom  row  of  plots)  during  high  solar  activity  2001 
Figure  4b  shows  the  same  but  during  low  solar  activity  2008.  The  modip  latitude  is  also  marked  at 
the  plots  by  white  line  as  the  low-latitude  plasma  bulk  follows  the  modip  frame.  At  both  solar  cycle 
conditions  it  is  seen  that  the  model  maps  reproduce  very  well  the  main  features  of  the  CODE  TEC 
maps;  some  quantitative  difference  is  evident  only  at  winter  model  map  where  the  equatorial 
anomaly  is  slightly  weaker  than  that  at  the  CODE  TEC  map  (Figure  5a,  upper  row  of  plots).  The 
hemispheric  asymmetry  of  the  equatorial  anomaly,  generated  mainly  by  the  thermospheric 
transequatorial  neutral  winds  blowing  from  the  summer  to  winter  hemisphere,  is  evident  well  at 
both  model  and  CODE  maps  but  only  at  low  solar  activity  (Figure  5b,  upper  row  of  plots).  The 
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hemispheric  symmetry  of  the  equatorial  anomaly  during  vernal  equinox  is  also  well  evident  at  both 
model  and  CODB  maps  during  high  and  low  solar  activity  (Figures  5a  and  5b,  bottom  rows  of  plots). 
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Figure  5a  Global  maps  in  geographical  frame  calculated  from  the  empirical  background  TEC  model  (left  column  of 
plots)  which  are  compared  with  the  CODE  TEC  maps  (right  column  of  plots)  at  12  UT  for  15  January  2001  (upper  row 
of  plots)  and  15  March  2001  (bottom  row  of  plots)  during  high  solar  activity  2001.  The  modip  latitude  is  also  marked 
by  white  line. 
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Figure  5b  The  same  as  in  Figure  5a,  but  during  low  solar  activity  2008. 


Figure  6a  presents  the  global  maps  from  the  background  TBC  model  (left  column  of  plots) 
compared  with  the  CODB  TBC  maps  (right  column  of  plots)  at  12UT  for  15  July  2001  (upper  row 
of  plots)  and  15  October  2001  (bottom  row  of  plots)  during  high  solar  activity  2001  while  Figure  6b 
shows  the  same  but  for  low  solar  activity  2008.  Again  the  degree  of  similarity  between  model  and 
CODB  TBC  maps  is  very  high.  Some  undervalue  of  the  model  TBC  is  seen  in  July  at  most  northern 
latitudes  at  both  high  and  low  solar  activity  (Figures  6a  and  6b,  upper  row  of  plots).  The  model 
however  reproduces  very  well  the  four  TBC  amplifications  seen  between  ~50°N  and  -dO^S  at  the 
most  western  longitudes  and  two  TBC  amplifications  at  most  eastern  longitudes  in  July  2008 
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(Figure  6b,  upper  row  of  plots).  It  is  worth  noting  that  both  model  and  CODB  maps  show  the 
following  features:  (i)  do  not  display  hemispheric  asymmetry  of  the  equatorial  anomaly  during  July 
at  both  solar  activity  conditions,  (ii)  the  January  TBC  is  larger  than  that  in  July  at  high  and  low  solar 
activity  (so  called  winter  anomaly),  and  (ii)  while  the  March  TBC  is  higher  than  the  October  one 
during  low  solar  activity  the  opposite  feature  is  evident  during  high  solar  activity.  We  remind  also 
that  all  the  above  mentioned  features  are  seen  at  12UT. 
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Figure  6a  Global  maps  in  geographical  frame  calculated  from  the  empirical  background  TEC  model  (left  column  of 
plots)  which  are  compared  with  the  CODE  TEC  maps  (right  column  of  plots)  at  12  UT  for  15  July  2001  (upper  row  of 
plots)  and  15  October  2001  (bottom  row  of  plots)  during  high  solar  activity  2001.  The  modip  latitude  is  also  marked  by 
white  line. 
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Figure  6b  The  same  as  in  Figure  5a,  but  during  low  solar  activity  2008. 


The  comparisons  for  the  middle  solar  activity  2004  are  presented  in  Figures  7a  and  7b;  in  this  case 
only  maps  for  July  at  different  UT  are  shown  in  order  to  trace  out  better  the  diurnal  variability  of 
the  low  latitude  plasma  bulk.  Figure  7a  shows  global  maps  from  the  background  TBC  model  (left 
column  of  plots)  which  are  compared  with  the  CODB  TBC  maps  (right  column  of  plots)  at  OOUT 
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(upper  row  of  plots)  and  06UT  (bottom  row  of  plots)  for  15  July  2004  while  Figure  7b  presents  the 
same  but  at  12UT  (upper  row  of  plots)  and  at  18UT  (bottom  row  of  plots).  Again  the  comparison 
shows  high  degree  of  similarity;  even  some  details  in  the  spatial  TBC  distribution  are  well 
reproduced  (see,  for  example.  Figure  7a). 


TEC  Model 
15  Jul  2004,  OOUT 


-180  -150  -120  -90  -60  -30  0  30  60  90  120  150  180 


Longitude  (deg) 


15  Jul  2004,  06UT 


-180  -150  -120  -90  -60  -30  0  30  60  90  120  150  180 


Longitude  (deg) 


TEC  Data 

15  Jul  2004,  OOUT 


-180  -150  -120  -90  -60  -30  0  30  60  90  120  150  180 


Longitude  (deg) 


15  Jul  2004,  06UT 


-180  -150  -120  -90  -60  -30  0  30  60  90  120  150  180 


Longitude  (deg) 


Figure  7a  Global  maps  in  geographical  frame  calculated  from  the  empirical  background  TEC  model  (left  column  of 
plots)  which  are  compared  with  the  CODE  TEC  maps  (right  column  of  plots)  at  00  UT  (upper  row  of  plots)  and  06  UT 
(bottom  row  of  plots)  for  15  July  2004  during  middle  solar  activity  (2004).  The  modip  latitude  is  also  marked  by  white 
line. 


15  July  2004,  12UT 
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15  July  2004,  18UT 
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Figure  7b  The  same  as  in  Figure  7a,  but  at  12  UT  (upper  row  of  plots)  and  at  18  UT  (bottom  row  of  plots). 

The  above  shown  comparisons  indicate  that  the  empirical  background  TEC  model  can  reproduce 
very  well  the  main  spatial-temporal  variability  of  the  CODE  TEC  maps.  Each  empirical  model 
needs  to  be  statistically  evaluated.  A  detailed  description  of  the  model  error  will  be  done  in  the  next 
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section.  Here  however  only  the  main  statistics  based  on  the  entire  data  set  will  be  presented.  It  has 
been  accepted  that  the  mean  (systematic)  error  {ME),  root  mean  squares  error  (RMSE)  and  the 
standard  deviation  error  (STDE)  are  the  main  error  characteristics  of  each  model.  They  are  defined 
as: 


ME=jZ(TEC,^-TEC^) 


emseJ^XItec,^  -tecJ' 


(4) 


STDE  ^^RMSE^-ME^ 

The  application  of  (4)  to  all  data  for  the  considered  period  of  time  (1  January  1999-31  December 
2011)  gives  the  following  errors:  ME=0.003  TBCU,  i.e.  the  model  has  practically  zero  systematic 
error.  In  this  case  RMSE=STDE=33S1  TBCU.  In  order  to  have  an  idea  if  such  errors  are  high  or 
low  we  compare  this  model  with  the  similar  to  some  extend  new  global  TBC  model  built  recently 
by  Jakowski  et  al.  (2011)  and  called  NTCM-GF  model.  The  statistical  assessments  of  the  NTCM- 
GL  model  are:  ME=-0.3  TBCU  and  RMSE=1.5  TBCU.  Hence,  the  errors  of  the  present  background 
TBC  model  are  significantly  smaller  than  those  of  the  NTCM-GL  model.  We  have  to  note  that 
nevertheless  that  both  models  are  climatological,  i.e.  they  describe  the  average  behavior  under  quiet 
geomagnetic  conditions,  the  TEC  model  constructed  by  Jakowski  et  al.  (2011)  needs  only  12 
coefficients  and  a  few  empirically  fixed  parameters  for  describing  a  broad  spectrum  of  TEC 
variation  at  all  levels  of  solar  activity.  We  however  do  not  consider  the  large  number  of  coefficients 
in  the  present  background  TEC  model,  4374,  as  its  weak  point.  They  are  calculated  only  once  and 
are  fixed  later  at  the  model  applications. 


1.4  Concluding  comments  and  on-line  implementation  of  the  background  TEC 
model 

A  global  background  TEC  model  is  built  on  the  basis  of  full  13  years  (1999-2011)  of  CODE  TEC 
data  (Mukhtarov  et  al.,  2013a).  The  model  describes  the  climatological,  i.e.  under  quiet 
geomagnetic  conditions,  behaviour  of  the  ionosphere  and  can  be  used  for  long-term  prediction.  For 
this  purpose  at  given  day  of  the  year,  geographic  location  and  UT  the  model  needs  as  input 
parameters  only  the  predicted  level  of  solar  activity  (F10.7  is  used  here  as  a  proxy  for  solar  activity). 
The  model  maps  are  arrayed  in  terms  of  the  coordinate  system  of  geographical  latitude  from  -90°  to 
90°  at  each  5°  and  longitude  from  -180°  to  180°  at  each  5°. 

The  model  describes  very  well  such  structures  as  the  WSA  (Figure  3)  and  the  well-known  WN4  and 
WN3  longitude  structures  (Figure  4).  This  was  possible  mainly  because  of  the  nonmigrating  tides 
and  SPW  inclusion  in  describing  the  diurnal  variability  of  the  TEC.  The  presented  comparison 
between  the  model  and  CODE  TEC  maps  at  different  solar  activities  and  seasons  (Figures  5,  6  and 
7)  also  demonstrated  high  degree  of  similarity.  The  zero  systematic  error  and  its  low  RMSE  (3.387 
TECU)  provides  the  model  significant  advantage  over  the  other  similar  models. 

The  present  background  model  can  be  used  for  both  science  and  applications.  In  science  the  model 
can  be  utilized  as  a  background  condition  on  the  basis  of  which  the  perturbations  can  be  estimated. 
It  is  particularly  useful  for  investigating  the  geomagnetic  perturbations,  or  ionospheric  disturbances 
related  to  the  sudden  stratospheric  warmings,  by  incoherent  scatter  radars  where  the  measurements 
are  available  only  for  several  days;  in  this  case  the  background  condition  described  by  the  monthly 
median  TEC  values  cannot  be  determined.  This  model  can  be  useful  for  numerous  single  frequency 
GPS  applications  which  need  additional  information  for  mitigating  the  ionospheric  propagation 
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error.  Such  GNSS  users  can  be  provided  by  adequate  ionospheric  corrections  obtained  by  this 
autonomous  ionospheric  TEC  model. 


This  model  will  be  used  for  long-term  prediction  of  the  TEC.  The  online  software  has  been  already 
implemented  on  the  website  (http://www.geophys.bas.bg/tec/page/tec_index_en.htm),  but  it  is  still 
in  testing  mode.  Figure  8  shows  its  block-chart.  It  will  be  officially  implemented  at  the  institute 
website  when  the  paper  Mukhtarov  et  al.  (2013a)  is  published  (it  has  been  accepted  for  publication 
in  J.  Geophys.  Res.  -  Space  Physics). 


Long  -  time  forecast 


Current  and  forecast 
Solar  activity 


Figure  8  The  block  chart  of  the  TEC  model  for  long-term  prediction 


2.  Statistical  evaluation  of  the  global  background  TEC  model 
2.1  General  Evaluation  of  the  Background  TEC  Model 

An  important  aspect  of  the  model  development  process  is  the  evaluation  of  model  performance 
comprehensively  and  objectively.  This  means  that  we  have  to  represent  an  objective  and  meaningful 
description  of  the  model's  ability  to  reproduce  reliable  observations  precisely  or  accurately,  i.e.  to 
determine  the  extent  to  which  model-predicted  events  approach  a  corresponding  set  of  reliable 
observations  To  gain  even  further  insight  into  the  nature  and  sources  of  the  model’s  error  variable 
we  study  in  detail  the  solar,  seasonal  and  diurnal  variability  of  the  error  and  on  the  base  of  the 
obtained  results  we  will  present  an  error  model  as  well. 
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In  Section  1  we  have  already  presented  the  overall  statistical  assessment  of  the  model  based  on  the 
entire  data  set.  The  model  performance  has  been  represented  by  the  mean  (systematic)  error  (MF), 
root  mean  squares  error  (RMSE)  and  the  standard  deviation  error  (STDE)  calculated  with  the 
expressions  (4)  in  the  Section  1 .  It  has  been  found  that  the  background  model  fits  to  the  CODE  TEC 
input  data  with  a  zero  systematic  error  and  a  RMSE=STDE=3. 3S7  TECU.  The  empirical  probability 
density  distribution  of  the  model’s  error  is  shown  in  Figure  9a  (black  line).  It  is  almost  a  symmetric 
function  and  bears  a  resemblance  in  some  way  to  the  Laplace  distribution,  shown  in  Figure  la  by 
red  line  (calculated  at  the  same  mathematical  expectation  and  variance  as  the  empirical  one),  but 
has  also  significant  differences  particularly  around  errors  close  to  zero.  The  confidence  limits  of  the 
error  at  a  given  probability  are  determined  empirically  by  numerical  integration  of  the  probability 
density  function  shown  by  black  line  in  Figure  9a.  Figure  9b  shows  the  probability  for  obtaining  a 
given  error  expressed  in  times  STDE  (i.e.  the  empirical  error  function).  Figure  9b  reveals  that  the 
90%  probability  corresponds  to  an  error  interval  from  -\.5STDE  to  \.5STDE,  i.e.  from  about  -5  to  5 
TECU.  This  means  that  there  is  a  90%  probability  that  deviations  larger  than  5  TECU  between  the 
model  and  the  CODE  TEC  data  would  not  occur. 


(a) 


-15  -12  -9  -6  -3  0  3  6  9  12  15 

Error  (TECU) 

(b) 


Figure  9  (a)  Empirical  probability  density  distribution  of  the  model’s  error  (black  line)  compared  with  the  respective 
(calculated  at  the  same  mathematical  expectation  and  variance  as  the  empirical  one)  Laplace  distribution  (red  line);  (b) 
Probability  for  obtaining  a  given  error  expressed  in  times  STDE. 

The  overall  statistics  of  the  model  error  can  be  defined  more  precisely  by  showing  its  dependence 
on  LT  and  modip  latitude.  Figure  10  shows  the  mean  (systematie)  error  {ME)  dependence  on  modip 
latitude  and  LT.  It  is  seen  that  it  reaches  the  largest  values  of  +0.7  TECU  (insignificant  error) 
mainly  at  low-  and  low-middle  latitudes.  The  ME  variability  reflects  the  faet  that  the  fifth  harmonics 
of  the  solar  day  (4.8-hour  tidal  component)  is  not  included  in  the  background  TEC  model. 
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Figure  10  Mean  (systematic)  error  dependence  on  modip  latitude  and  LT 


Figure  11  Dependence  of  RMSE  (upper  plot;  the  contour  distance  is  0.5  TECU)  and  relative  RMSE  (bottom  plot,  the 
contour  distance  is  0.025)  on  modip  latitude  and  LT. 

Figure  1 1  shows  the  RMSE  distribution  (upper  plot;  the  contour  distance  is  0.5  TBCU)  with  respect 
to  the  modip  latitude  and  FT.  The  largest  errors  are  obtained  around  sunrise  (~8  FT)  and  sunset 
(~18  LT).  While  the  sunrise  errors  maximise  above  the  equator  (and  this  is  normal  because  of  the 
absence  of  equatorial  ionization  anomaly  those  at  sunset  maximise  not  only  above  the  equator  but 
also  at  +30°  modip  latitude,  i.e.  at  the  equatorial  anomaly  crests.  The  errors  at  the  northern 
equatorial  anomaly  crest  are  slightly  larger  than  those  at  the  southern  crest.  This  result  could  be 
connected  with  the  asymmetric  behaviour  of  the  migrating  diurnal  (DWl)  and  semidiurnal  (SW2) 
tidal  components  seen  at  the  left  column  of  plots  in  Figure  2  of  Section  1.  Both  tidal  components 
are  stronger  in  the  Northern  Hemisphere  (NH)  however  this  asymmetry  for  the  DWl  is  better 
expressed  at  high  solar  activity  while  for  the  SW2  is  well  seen  at  low  solar  activity.  This 
asymmetric  tidal  behaviour  is  not  well  described  by  the  background  TBC  model.  The  bottom  plot  of 
Figure  1 1  presents  the  relative  RMSE  distribution  (the  contour  distance  is  0.025)  with  respect  to  the 
modip  latitude  and  LT.  It  is  seen  that  the  relative  error  is  small  during  the  daytime  (between  8  and 
18  LT)  everywhere;  it  is  particularly  small,  not  more  than  5%,  between  12  and  16  LT  around  the  dip 
equator.  The  largest  relative  errors  of  -30%  can  be  distinguished  between  2  and  4  LT  above  the  dip 
equator  and  between  18  and  6  LT  (night-time)  at  around  -60°  modip  latitude.  The  former  largest 
relative  error  area  is  most  probably  related  to  the  temperature  and  broad  plasma  anomalies  {Huang 
at  ai,  2012)  observed  above  the  equator  around  and  after  midnight;  they  are  considered  as  part  of 
the  tidal  pattern.  The  latter  largest  relative  error  area  is  surely  related  to  the  WSA  discussed  in  the 
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Section  1.  This  means  that  nevertheless  that  the  background  TBC  model  is  able  to  model  the  WSA 
its  variability  is  a  source  of  large  errors.  Some  increase  of  the  relative  error  is  seen  also  at  about 
+30°  modip  latitude  between  20  and  23  LT  which  is  most  probably  related  to  oscillations  in  the 
equatorial  evening  pre-reversal  electric  field  (F-region  vertical  drift)  and  their  effect  on  the 
variability  of  the  plasma  irregularities. 


2.2  Basic  Approach  of  the  Error  Model  Construction 

In  order  to  asses  the  dependence  of  the  error  on  the  solar  activity,  seasons  and  LT  we  have  to 
demonstrate  how  the  model’s  error  changes  at  different  conditions.  For  this  purpose  we  calculated 
the  monthly  mean  values  of  the  RMSE  for  the  considered  period  of  time,  1  January  1999  -  31 
December  2011.  The  left  column  of  plots  in  Figure  12  shows  the  modip  latitude-time  cross  sections 
of  the  monthly  mean  RMSE  at  different  LT:  OOLT  (most  upper  plot),  08LT  (second  from  above 
plot),  12LT  (third  from  above  plot)  and  18LT  (bottom  plot).  It  is  clearly  evident  that  the  model’s 
errors  are  larger  during  high  solar  activity  (i.e.  solar  cycle  dependence),  at  equinoxes  (i.e.  have 
seasonal  dependence)  and  they  depend  on  LT.  All  the  above  mentioned  dependences  are  very 
similar  to  those  of  the  TBC  itself.  Similarly  to  the  TBC  model  here  also  the  time  scales  of  the  error 
variability  related  to  the  solar  cycle,  season  and  LT  are  very  different.  Therefore  for  building  the 
error  model  we  use  the  same  approach  as  that  applied  in  constructing  the  background  TBC  model 
and  the  RMSB  can  be  represented  as: 

RMSE  {E,Kp,  month,  (F,  )'R2  {month)%  (LF)  (5) 


Similarly  to  the  TBC  model  here  also  the  above  right  hand  side  unknown  functions  Wk  (k=l,2,3)  can 
be  represented  by  their  series  expansions;  Wj  can  be  expanded  in  Taylor  series,  while  W2  and  W 3, 
which  are  periodic  functions  with  periods  respectively  a  year  (12  months)  and  a  solar  day  (24  hours), 
can  be  expanded  in  Fourier  series.  Therefore,  at  each  modip  latitude  the  error  model  can  be 
described  by  the  following  functions: 


RMSE{E,  Kp, month, LT)  =  +  a^E  +  a-^Kp  +  a^E^  +  a^EKp  +  u^Kf^) 
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The  expression  in  the  first  right  hand  bracket,  i.e.  the  Taylor  series  expansion  up  to  degree  of  2, 
represents  the  solar  activity  term  which  modulates  the  seasonal  and  diurnal  behavior  of  the  RMSE. 
Similarly  to  the  Section  1  here  also  F  is  the  solar  radio  flux  at  10.7  cm  wavelength  (F10.7)  and  Kf 
describes  the  linear  rate  of  change  of  F10.7.  The  seasonal  term  (expression  in  the  second  right  hand 
bracket)  includes  besides  the  yearly  mean  (bo)  also  4  subharmonics  of  the  year,  i.e.  annual, 
semiannual,  4-  and  3-month  components;  it  modulates  the  diurnal  behavior  of  the  RMSE.  In  this 
case  the  diurnal  variability  of  the  RMSE  model  (expression  in  the  third  right  hand  bracket)  is 
composed  however  only  by  two  terms:  daily  mean  RMSE  (cq)  and  a  term  describing  the  migrating 
tides.  This  is  due  to  the  fact  that  the  RMSE  depends  mainly  on  the  LT.  The  contribution  of  the 
migrating  tides  in  (6)  includes  4  subharmonics  of  the  solar  day,  i.e.  24-,  12-,  8-  and  6-hour 
components. 

The  error  model  described  by  (6)  contains  486  constants  (we  remind  that  (6)  is  applied  at  each 
modip  latitude)  and  they  are  determined  by  least  squares  fitting  techniques.  Similarly  to  the  TBC 
model  here  also  the  numbers  of  the  included  components  in  the  Taylor  and  Fourier  expansion  series 
are  defined  experimentally.  We  accepted  only  the  above  mentioned  solar,  seasonal  and  diurnal 
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components  because  the  addition  of  more  components  does  not  improve  significantly  the  error  of 
the  error  model. 


RMSE  tOR  LT 


12  24  36  48  60  72  84  96  108  120  132  144  156 


Model  RMSRmOT  X 


Month  Number  (start  January  1999)  Month  Number  (start  January  1999) 


Figure  12  (left  column  of  plots)  Modip  latitude-time  (months)  cross-sections  of  the  calculated  monthly  mean  RMSE  at  different  local 
times:  00  LT  (most  upper  plot),  08  LT  (second  from  above  plot),  12  LT  (third  from  above  plot)  and  18  LT  (bottom  plot);  (right 
column  of  plots)  the  same  as  the  left  column  of  plots  but  for  the  model  RMSE. 


The  error  model  offers  a  prediction  approach  on  the  basis  of  which  we  can  predict  the  RMSE 
depending  on  the  solar  activity,  season  and  LT.  Therefore,  we’ve  built  not  only  a  global  empirical 
background  TEC  model  but  also  a  global  prediction  of  the  model’s  error  at  different  solar,  seasonal 
and  LT  conditions.  The  overall  standard  deviation  of  the  predicted  RMSE  with  respect  to  the 
empirical  obtained  one  is  0.7  TECU.  The  right  column  of  plots  in  Figure  12  presents  the  same 
results  as  the  left  column  of  plots  but  for  the  model  RMSE.  The  detailed  comparison  between  the 
real  and  model  RMSE  reveals  some  important  features.  The  error  model  describes  very  well  the  real 
RMSE  at  00  and  08  LT;  it  is  able  to  reproduce  not  only  qualitatively,  but  also  quantitatively  solar 
and  seasonal  dependences  of  the  RMSE.  The  error  model  is  able  to  reproduce  even  the  hemispheric 
asymmetry  of  the  RMSE  well  seen  particularly  at  high  solar  activity;  it  is  larger  in  the  SH  at  00  LT 
and  in  the  NH  at  08  LT.  The  error  model  performance  at  12  and  18  LT  is  not  as  good  as  that  at  00 
and  08  LT.  The  model  has  not  been  able  to  reproduce  well  particularly  the  large  errors  seen  in  1999, 
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2001  and  2011.  Generally  however  the  error  model  describes  correctly  the  solar  and  seasonal 
dependences  of  the  RMSB  and  its  global  distribution. 

2.3  Application  of  the  Error  Model 

The  global  empirical  background  TBC  model,  described  in  Section  1,  offers  TEC  maps  which 
depend  on  geographic  coordinates  (5°x5°  in  latitude  and  longitude)  and  UT  at  given  solar  activity 
and  day  of  the  year.  The  error  model  however  does  not  depend  on  the  geographic  longitude  because 
only  the  contribution  of  the  migrating  tidal  components  is  considered  in  the  model.  In  this  way  the 
error  maps  depending  on  the  geographic  latitude  and  FT  at  given  solar  activity  and  month  of  the 
year  have  to  be  constructed.  The  conversion  of  the  modip  latitude  to  geographic  one  is  done  at  the 
Greenwich  meridian.  The  error  values  assigned  to  both  poles  are  the  same  as  those  at  the  highest 
northern  and  southern  latitudes. 

In  order  to  demonstrate  the  ability  of  the  error  model  to  reproduce  the  spatial-temporal  features  of 
the  real  RMSE  at  different  solar  activity  and  seasons  we  use  the  examples  given  in  the  Section  1 
(Figures  5,  6  and  7),  i.e.  we  will  compare  the  real  and  model  RMSE  for  2001,  2004  and  2008  as 
years  representing  high,  middle  and  low  solar  activity  and  months  January,  March,  July  and 
October  as  typical  winter,  spring,  summer  and  autumn  months. 
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Figure  13  Latitude-LT  cross  sections  of  the  real  (left  column  of  plots)  and  model  (right  column  of  plots)  RMSE  for  15 
January  (upper  row  of  plots),  15  March  (second  from  above  row  of  plots),  15  July  (third  from  above  row  of  plots)  and 
15  October  (bottom  row  of  plots)  at  high  solar  activity  2001. 
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Figure  13  presents  latitude-LT  cross  sections  of  the  real  (left  column  of  plots)  and  model  (right 
column  of  plots)  RMSE  for  15  January  (upper  row  of  plots),  15  March  (second  from  above  row  of 
plots),  15  July  (third  from  above  row  of  plots)  and  15  October  (bottom  row  of  plots)  at  high  solar 
activity  2001.  The  latitude-LT  distributions  of  the  real  and  model  RMSE  in  January  (upper  row  of 
plots)  are  very  similar  not  only  qualitatively  but  quantitatively  as  well.  As  usually  the  largest  errors 
are  seen  at  both  plots  around  sunrise  (~6-8  LT)  and  sunset  (-16-20  LT);  large  errors  are  found  at 
both  plots  also  near  20°N  mainly  during  the  daytime  (-6-20  LT)  and  above  the  equator  at  sunset. 
The  degree  of  similarity  between  the  real  and  model  RMSE  in  March  (second  from  above  row  of 
plots)  is  also  very  high;  at  both  plots  the  errors  are  symmetrically  distributed  with  respect  to  latitude 
of  -10°N  (because  of  difference  between  modip  and  geographic  latitudes).  Again  the  largest  errors 
at  both  plots  are  seen  near  sunrise  (-7  LT)  and  sunset  (-18  LT)  but  while  the  maximum  real  RMSE 
is  1 1  TBCU  that  of  the  model  RMSE  is  slightly  weaker,  i.e.  it  is  10  TBCU.  The  comparison  between 
summer  (July)  real  and  model  RMSE  (third  from  above  row  of  plots)  again  demonstrates  high 
degree  of  similarity,  qualitatively  and  quantitatively.  In  this  case  the  largest  errors  are  seen  only  in 
the  morning  hours  (-8-10  LT)  and  at  -20°N,  but  not  during  the  sunset.  The  qualitative  similarity 
between  the  real  and  model  RMSE  in  October  is  very  good  however  the  model  errors  are  smaller 
than  those  of  the  real  RMSE,  i.e.  they  are  8.1  TBCU  and  1 1  TBCU  respectively. 


Figure  14  The  same  as  Figure  5  but  at  low  solar  activity  2008. 
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Figure  14  shows  the  same  comparison  as  that  in  Figure  13  but  for  low  solar  activity  2008.  In  this 
case  both  real  and  model  RMSE  drastically  decrease.  The  model  describes  qualitatively  very  well 
the  latitude-LT  distribution  of  the  real  RMSE  at  all  months;  there  is  some  quantitative  difference 
mainly  during  the  equinoxes.  As  it  has  been  expected  the  largest  errors  are  found  at  equinoxes  both 
in  real  and  model  RMSE  but  the  error  model  underestimates  the  RMSE;  in  March  the  model  and  real 
RMSE  are  respectively  4  and  6  TBCU,  while  for  October  the  difference  is  smaller  and  they  are  3.2 
and  4.5  TECU.  Some  hemispheric  asymmetry  of  both  real  and  model  RMSE  is  seen  in  winter 
(January)  and  summer  (July)  here  as  well  but  it  is  not  predominantly  in  the  NH  as  it  was  at  high 
solar  activity  2001  (Figure  13).  In  January  both  real  and  model  RMSE  at  sunrise  and  morning  hours 
are  stronger  in  the  NH  while  at  afternoon  and  sunset  hours  they  are  larger  in  the  SH.  The  opposite 
asymmetry  is  seen  in  July. 


Figure  15  The  same  as  Figure  5  but  at  middle  solar  activity  2004. 


Figure  15  presents  a  comparison  between  real  and  model  RMSE  maps  at  middle  solar  activity  2004. 
It  is  seen  that  both  the  real  and  model  RMSE  at  all  months  are  between  those  at  high  (Figure  13)  and 
low  (Figure  14)  solar  activity,  as  it  is  expected  in  advance.  At  all  months  the  real  and  model  largest 
values  are  similar  but  in  March  they  are  almost  the  same,  i.e.  6.9  and  6.8  TECU.  The  largest 
difference  is  seen  in  July  when  the  maximum  real  RMSE  is  4.5  TECU  while  the  model  one  is  3.6 
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TBCU.  During  the  daytime  almost  at  all  months  the  NH  errors  are  larger  than  those  in  the  SH;  only 
in  January  both  the  real  and  model  RMSE  distribution  and  the  real  RMSE  in  March  are  more 
hemispheric  symmetric.  Similarly  to  high  solar  activity  here  also  the  increase  of  RMSE  during 
daytime  and  at  ~20'’N  is  a  consequence  of  the  hemispheric  asymmetry  of  the  diurnal  components 
DWl  and  SW2  contribution  to  the  equatorial  anomaly  which  is  not  well  reproduced  by  the 
background  TBC  model.  While  the  night  time  (~2-4  FT)  amplification  of  the  RMSE  in  March  is 
comparatively  well  reproduced  that  observed  in  January  is  underestimated. 

2.4  Concluding  comments 

A  detailed  statistical  evaluation  of  the  global  empirical  background  TBC  model,  presented  in 
Section  1,  is  done  in  Section  2.  The  model  performance  has  been  described  by  its  mean  (systematic) 
error  {ME),  root  mean  squares  error  {RMSE)  and  the  standard  deviation  error  {STDE).  It  was  found 
that  the  background  model  fits  to  the  CODB  TBC  input  data  with  a  zero  systematic  error  and  a 
RMSE=STDE=3.3S1  TBCU.  Based  upon  this  overall  error  measures  we  may  confidently  conclude 
that  this  model  is  able  to  reproduce  accurately  the  CODE  TEC  input  data.  It  was  found  that  the 
empirical  probability  density  distribution  (Figure  9a)  resembles  more  the  Laplace  than  the  normal, 
or  Gaussian,  distribution.  This  result  could  be  probably  related  to  the  non-Gaussian  statistics  of  the 
ionospheric  irregularities.  The  empirical  error  function  shown  in  Figure  9b  revealed  that  there  is 
only  10%  probability  that  deviations  larger  than  5  TBCU  between  the  model  and  the  CODE  TEC 
data  would  occur.  The  modip  latitude-LT  distributions  of  the  model’s  error  showed  predominantly 
known  features,  as:  (i)  the  small  ME  observed  mainly  at  low  latitudes  reflects  the  fact  that  the  fifth 
harmonics  of  the  solar  day  (4.8-hour  tidal  component)  is  not  included  in  the  background  TEC  model 
(Figure  10);  (ii)  the  RMSE  are  large  at  sunrise  and  sunset  time  (Figure  11  upper  plot),  and  (iii)  the 
relative  RMSE  amplifications  shown  in  the  bottom  plot  of  Figure  1 1  are  related  to  comparatively 
stable  ionospheric  anomalies  which  are  present  at  some  local  times  and  latitudes  (as  WSA,  broad 
plasma  anomaly  after  midnight  and  evening  pre-reversal  plasma  irregularities  at  equatorial 
latitudes). 

To  gain  further  insight  into  the  nature  and  sources  of  the  model’s  error  variable  we  studied  in  detail 
the  solar,  seasonal  and  diurnal  variability  (LT)  of  the  model’s  error.  On  the  base  of  the  obtained 
results  we  built  an  error  model  (Mukhtarov  et  al.,  2013b).  It  actually  offers  a  prediction  approach  on 
the  basis  of  which  we  can  predict  the  RMSE  depending  on  the  solar  activity,  season  and  LT.  The 
error  model  was  constructed  by  using  a  similar  approach  to  that  of  the  background  TEC  model  itself. 
Similarly  to  the  TEC  model  here  also  the  time  scales  of  the  error  variability  related  to  the  solar 
cycle,  season  and  LT  are  very  different.  Then  the  spatial-temporal  variability  of  the  RMSE  was 
presented  as  a  multiplication  of  three  separable  functions  (as  it  is  shown  in  (6)).  The  solar  cycle  and 
seasonal  dependences  of  the  RMSE  are  described  in  the  same  way  as  in  the  background  TEC  model. 
The  Taylor  series  expansion  up  to  degree  of  2  represents  the  solar  activity  function  while  the 
seasonal  function  includes  the  contribution  of  4  subharmonics  of  the  year,  i.e.  annual,  semiannual, 
4-  and  3-month  components.  The  RMSE  depends  mainly  on  the  LT  and  due  to  this  its  diurnal 
variability  is  described  only  by  the  migrating  tides;  four  subharmonics  of  a  solar  day,  24-,  12-,  8- 
and  6-hour  components  are  included  in  the  error  model.  It  contains  486  constants  which  have  been 
determined  by  least  squares  fitting  techniques.  The  overall  standard  deviation  of  the  predicted 
RMSE  with  respect  to  the  empirical  one  is  0.7  TECU.  The  detailed  comparisons  between  real  and 
model  RMSE  shown  in  Figures  13,  14  and  15  clearly  demonstrate  that  the  error  model  describes 
correctly  and  precisely  the  spatial-temporal  variability  of  the  RMSE. 

It  is  important  to  note  that  these  two  sections  present  not  only  a  global  empirical  background  TEC 
model  (Mukhtarov  et  al.,  2013a)  but  also  a  global  prediction  of  the  model’s  error  at  different  solar, 
seasonal  and  LT  conditions  (Mukhtarov  et  al.,  2013b).  At  given  solar  activity  and  day  of  the  year 
the  background  TEC  model  offers  TEC  maps  which  depend  on  geographic  coordinates  (5°x5°  in 
latitude  and  longitude)  and  UT.  The  error  model  offers  a  prediction  approach  on  the  basis  of  which 
the  RMSE  depending  on  the  solar  activity,  season  and  LT  can  be  predicted. 
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3.  Development  of  global  empirical  model  of  TEC  response  to  geomagnetic 
activity 

The  model  is  constructed  on  the  basis  of  the  same  CODB  TBC  data  set.  The  geomagnetic  activity  is 
defined  by  the  global  .^f^-index.  The  .^f^-index  data  are  downloaded  from  the  Space  Physics 
Interactive  Data  Resource  (SPIDR),  Boulder,  Colorado  for  the  considered  period  of  time.  The  Kp 
value  at  every  hour  is  used  in  this  study  as  it  is  obtained  by  linear  interpolation  of  the  three-hour  Kp 
values. 

In  the  model  the  TBC  response  to  the  geomagnetic  activity  is  investigated  by  the  relative  deviation 
of  the  TBC  defined  as:  rTBC  =  (TBCobs-TBCmed)/TBCmed-  The  terms  TBCobs  and  TBCmed  represent 
observed  and  median  TBC  values  respectively  at  a  given  hour.  In  the  present  study  we  use  sliding 
medians  defined  by  a  15-day  moving  window  and  the  median  value  is  assigned  to  the  last  day  of  the 
window,  i.e.  to  the  15*  day  of  the  window.  We  use  such  type  of  “one-sided”  median  approach 
because  in  this  case  the  model  is  built  in  a  way  as  it  will  be  used  for  short-term  prediction  (usually 
24  hours  ahead).  In  other  words  we  try  to  predict  the  correction  to  the  15-day  median  values  for 
each  hour  of  the  prediction  period.  The  window  length  of  15  days  is  used  because  the  preliminary 
examinations  showed  both  strong  depression  of  the  ~27-day  rTEC  variability  due  to  analogous 
variations  in  the  EUV  radiation  and  an  insignificant  effect  on  the  9-  and  13.5-day  recurrent 
geomagnetic  activity  particularly  strong  in  declining  phase  of  the  solar  activity.  By  considering  the 
characteristic  rTEC  the  effect  of  the  regular  seasonal,  diurnal  and  solar  changes  is  removed  from  the 
TEC  variability.  The  data  are  grouped  into  12-month  bins  as  every  bin  contains  all  the  available 
hourly  data  within  the  respective  month  of  the  year. 

3.1  Cross-correlation  analysis  between  rTEC  and  Kp-index  and  its  theoretical 
substantiation 

The  investigations  on  the  foF2  response  to  the  geomagnetic  activity,  presented  by  Muhtarov  and 
Kutiev  (1998)  and  Kutiev  and  Muhtarov  (2003),  indicated  that  this  is  a  delayed  response.  The 
authors  expressed  the  delay  in  terms  of  the  time  constant  of  their  cross-correlation  function  and 
found  a  time  delay  constant  of  about  18  hours.  Then  the  first  step  in  this  study  is  to  calculate  the  2D 
cross-correlation  functions  between  the  K^-index  and  rTEC. 

3.1.1  Empirical  cross-correlation  functions 

The  effect  of  geomagnetic  activity,  described  by  the  K^-index,  on  the  rTEC  variability  is 
investigated  by  2D  cross-correlation  analysis  between  both  parameters.  In  the  case  of  building 
global  TEC  model  we  expect  that  the  cross-correlation  function  will  depend  on  the  season,  modip 
latitude,  longitude  and  LT.  Due  to  these  dependences  three  different  types  of  2D  cross-correlation 
functions  are  calculated:  (i)  longitude-time  lag;  (ii)  LT-time  lag,  and  (iii)  modip  latitude-time  lag. 
The  2D  cross-correlation  functions  are  calculated  for  each  month  of  the  year  because  they  depend 
on  the  season  as  well.  Only  some  examples  of  the  above  mentioned  three  types  of  the  2D  cross¬ 
correlation  functions  will  be  shown  here  through  which  the  main  features  of  the  geomagnetic  effects 
on  the  rTEC  can  be  demonstrated. 

Figure  16  presents  2D  (longitude-time  lag)  cross-correlation  functions  calculated  between  the  rTEC 
and  Kp-index  for  January  at  different  modip  latitudes:  equator  (upper  most  plot),  +20°  (upper  row  of 
plots),  +40°  (middle  row  of  plots)  and  +70°  (bottom  row  of  plots).  The  2D  cross-correlation 
functions  from  both  hemispheres  are  shown  in  order  to  demonstrate  the  seasonal  dependence  of  the 
TEC  response  to  the  geomagnetic  activity;  winter  in  the  Northern  Hemisphere  (NH)  and  summer  in 
the  Southern  Hemisphere  (SH).  The  time  lag  up  to  72  hours  is  shown  because  in  general  the 
response  is  composed  by  two  phases,  positive  and  negative  with  different  duration  and  time  delay. 
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Figure  16  Two-dimensional  (longitude-time  lag)  cross-correlation  functions  calculated  between  the  rTEC  and  ^Tp-index 
for  January  at  equator  (upper  most  plot),  +20°  (upper  row  of  plots),  +40°  (middle  row  of  plots)  and  +70°  (bottom  row  of 
plots);  the  zero  time  lag  is  marked  by  tick  black  line. 

Some  main  features  of  the  cross-correlation  can  be  distinguished  from  Figure  16:  (i)  the  rTEC 
response  to  the  Kp-index  shows  clear  longitude  and  even  some  wave-like  dependence;  in  the  NH 
(winter)  a  wavenumber  3  response  can  be  seen  while  in  the  SH  (summer)  and  over  the  equator  in 
general  wavenumber  2  can  be  clarified;  (ii)  two  types  of  response,  positive  and  negative,  can  be 
seen  at  all  plots  (the  zero  time  lag  is  marked  by  tick  black  line);  first  the  response  is  positive  at  all 
modip  latitudes  except  that  at  70°S  (i.e.  summer  high  latitudes  where  the  cross-correlation  reaches 
maximum  of  -0.5  with  an  average  time  lag  of  6  hours),  and  then  it  is  negative;  above  the  equator  it 
is  mainly  positive;  (iii)  the  maximum  positive  cross-correlation  of  -1-0.4  is  seen  at  winter  high 
latitudes  (70'’N)  which  is  reached  in  the  frame  of  1-3  hours  after  the  maximum  Kp-index;  with 
decreasing  the  modip  latitude  the  time  lag  for  reaching  maximum  increases  and  above  the  equator  it 
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is  on  the  average  after  ~9  hours.  The  results  from  Figure  16  reveal  that  in  the  global  rTBC  model  a 
dependence  on  the  longitude  has  to  be  included. 


Figure  17  Two-dimensional  (LT-time  lag)  cross-correlation  functions  calculated  between  the  rTEC  and  ^^-index  for 
September  at  equator  (upper  most  plot),  +20°  (upper  row  of  plots),  +40°  (middle  row  of  plots)  and  +70°  (bottom  row  of 
plots);  the  zero  time  lag  is  shown  by  tick  black  line. 

Figure  17  shows  2D  (FT -time  lag)  cross-correlation  functions  calculated  between  the  rTEC  and  Kp- 
index  for  September  at  the  same  modip  latitudes  as  in  Figure  16.  In  this  case  an  equinoctial  month 
is  shown,  autumnal  month  in  the  NH  and  vernal  month  in  the  SH.  Again  two  types,  positive  and 
negative,  of  the  rTEC  response  are  seen  at  all  plots.  The  following  cross-correlation  features  can  be 
summarized  from  Figure  2:  (i)  low-latitude  rTEC  response  is  mainly  positive;  the  negative  response 
is  reached  at  large  time  lags;  the  maximum  positive  correlations  are  obtained  between  8-10  LT  and 
18-20  LT  with  an  average  time  lag  of  6-9  hours  for  the  tropics  and  9-12  hours  above  the  equator; 
(ii)  middle  (+40°)  latitude  rTEC  response  clearly  indicates  first  positive  and  then  negative  phases 
with  different  durations;  the  maximum  positive  coefficients  are  reached  during  the  day-time  10-12 
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LT  and  around  midnight  with  time  lags  of  ~3  hours  for  the  NH  and  between  3  and  6  hours  for  the 
SH;  the  negative  response  is  stronger  for  the  NH  than  that  in  the  SH  reaching  maximum  around  2-4 
LT  and  ~18  LT;  (iii)  high  (+70°)  latitude  rTBC  response  is  defined  by  negative  phase  during  the  day 
(6-20  LT)  and  a  positive,  almost  instantaneous  (zero  time  lag),  response  during  the  night;  the 
negative  response  in  the  NH  is  stronger  than  that  in  the  SH.  The  results  from  Figure  17  reveal  that 
in  the  global  rTEC  model  a  dependence  on  the  LT  has  to  be  included  as  well. 
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Figure  18  Two-dimensional  (modip  latitude-time  lag)  cross-correlation  functions  calculated  between  the  rTEC  and  Kp- 
index  for  months:  January  (upper  row  of  plots),  March  (second  row  of  plots),  June  (third  row  of  plots)  and  September 
(bottom  row  of  plots)  and  at  0°E  (left  column  of  plots)  and  270°E  (right  column  of  plots);  the  zero  time  lag  is  shown  by 
tick  black  line. 


Figure  18  presents  the  2D  (modip  latitude-time  lag)  cross-correlation  functions  calculated  between 
the  rTEC  and  ^^-index  for  different  months:  January  (upper  row  of  plots),  March  (second  row  of 
plots),  June  (third  row  of  plots)  and  September  (bottom  row  of  plots)  and  at  two  longitudes:  0°E 
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(left  column  of  plots)  and  270°B  (right  column  of  plots).  Again  two  types  of  rTBC  response  can  be 
seen  with  different  duration  and  time  lag  which  depends  on  the  season  and  modip  latitudes.  In 
general:  (i)  the  tropical  latitudes  at  all  seasons  have  positive  response  with  large  time  lags;  (ii)  while 
the  winter  high  latitude  rTBC  has  first  positive  response  with  short  time  lags  and  then  weak 
negative  one  the  summer  rTBC  demonstrates  only  negative  response;  (iii)  the  middle  (up  to  +50°) 
latitude  rTBC  response  shows  first  a  weak  positive  response  with  short  time  lags  and  then  stronger 
negative  response  with  large  time  lags.  Considering  all  months  of  the  rTBC  responses  the  following 
feature  can  be  drawn:  the  rTEC  response  in  March/April  is  close  to  the  winter  response  in  the 
NH/SH  and  to  the  summer  one  in  the  SH/NH  while  the  rTBC  response  in  September/October  is 
close  to  the  summer  response  in  the  NH/SH  and  to  the  winter  response  in  the  SH/NH.  Similar  result 
but  only  for  the  NH  was  found  in  the  regional  TEC  model  (Andonov  et  al.,  201 1). 


3.1.2  Theoretical  cross-correlation  function  between  rTEC  and  Kp-index 

The  main  conclusion  from  all  types  of  cross-correlation  functions,  shown  in  Figures  16,  17  and  18, 
is  the  existence  of  two  types  of  the  rTEC  response,  positive  and  negative,  with  different  durations 
and  time  lags.  Both  responses  depend  on  the  longitude,  modip  latitude,  season  and  LT.  The  cross¬ 
correlation  results  can  be  used  for  supporting  the  use  of  two  different  time  constants  in  building  the 
global  empirical  rTEC  model  in  a  way  as  it  has  been  already  done  in  the  regional  rTEC  model  over 
American  sector  (Andonov  et  al.,  2011).  The  use  of  two  different  time  constants  hints  for  the 
simultaneous  action  of  at  least  two  different  processes  that  define  the  rTEC  response  to 
geomagnetic  activity.  The  existence  of  at  least  two  processes  is  considered  in  Mukhtarov  and 
Pancheva  (2012)  where  the  ionospheric  response  to  the  high  speed  solar  wind  streams  is  studied  by 
using  the  COSMIC  electron  density  measurements. 

A  method  for  modeling  the  cross-correlation  function  between  the  relative  foF2  and  the 
geomagnetic  index  is  described  in  Muhtarov  et  al.  (2002)  where  the  delayed  response  is  represented 
by  a  linear  filter  from  the  first  order.  Similar  approach  we  apply  in  this  study  in  order  theoretically 
to  base  the  use  of  two  different  time  constants  in  establishing  the  global  rTEC  model  response  to  the 
geomagnetic  activity. 

If  we  assume  that  the  temporal  variability  of  the  geomagnetic  index  can  be  described  as  a  stationary 
random  process  (for  simplicity  noted  as  x(t))  while  the  rTEC  (noted  here  as  y(t))  is  a  result  of  the 
converting  the  geomagnetic  activity  by  two  independent  linear  filters  from  the  first  order,  then: 

y{t)^a°^h,  {4)x{t  -  ^)d^  +  U) 


The  transition  functions  of  the  both  filters  can  be  denoted  as: 
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Thus  the  ionospheric  response  is  represented  as  a  linear  combination  of  two  delayed  responses  with 
time  delay  constants  respectively  T /  and  T2  and  coefficients  of  proportionality  respectively  ct  and  (3. 
Actually,  the  ionospheric  response  is  not  a  deterministic  process  hence  the  above  mentioned 
quantities  have  to  be  interpreted  as  characteristics  of  the  most  probable  response  at  given  conditions. 

The  cross-correlation  function  between  the  processes  v  and  y  is  described  by  the  interrelations  of 
Wiener-  Lee: 
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(r)  =  j  {ah,  {^) + {^))r^^  {t  -  4)d^ 


(9) 


The  autocorrelation  function  of  the  geomagnetic  activity  can  be  represented  with  sufficient  accuracy 
[Muhtarov  et  al,  2002]  by  an  exponential  function: 
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where  the  magnitude  of  the  logarithmic  decrement  of  Rxx(t)  is  approximately  14  h  (Muhtarov  et  al., 
2002).  Having  in  mind  the  above  mentioned  assumptions  then  the  cross-correlation  function  can  be 
expressed  as: 
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At  x=0  both  expressions  and  their  first  derivatives  become  equal.  Figure  19  shows  a  comparison 
between  the  theoretical  (described  by  formula  (11))  cross-correlation  function  calculated  for  Ti=12 
h,  72=32  h,  CIJ=1,  /}=-!  and  for  simplicity  the  variance  of  the  geomagnetic  activity  is  accepted  to  be 
1  (upper  plot)  and  the  empirical  cross-correlation  function  between  Kp-index  and  rTEC  for  August 


0.3 
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Figure  19  Comparison  between  the  theoretical  cross-correlation  function  (upper  plot)  calculated  for  77=12  h,  72=32  h, 
0=1,  j3=-l  and  for  simplicity  the  variance  is  accepted  to  be  1,  and  the  empirical  cross-correlation  function  between  the 
i^^-index  and  rTEC  for  August  and  at  (40°N,  0°E)  (bottom  plot) 
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and  at  (40°N,  0°B)  (bottom  plot).  It  is  seen  that  the  main  features  of  the  two  cross-correlation 
functions  are  very  similar:  a  positive  response  with  small  time  lag  followed  by  a  negative  response 
with  longer  time  lag.  The  assumption  that  the  sum  response  is  shaped  by  two  responses:  a  positive 
with  small  time  constant  and  negative  one  with  three  times  longer  time  constant  is  set  in  the  model. 
In  this  way  at  the  range  of  positive  time  lag  a  near  area  of  positive  correlation  and  a  distant  area  of 
negative  correlation  are  formed.  In  the  presence  of  only  one  process  it  is  impossible  both  positive 
and  negative  correlations  to  be  obtained.  The  investigation  of  the  relative  foF2  response  to  the 
geomagnetic  activity  in  summer  at  the  middle  latitudes  (Muhtarov  et  al.,  2002;  Kutiev  and 
Muhtarov,  2003)  indicated  that  the  response  is  negative.  In  this  study  however  the  rTBC  response  at 
the  same  conditions  is  composed  by  positive  and  negative  responses.  This  means  that  the  positive 
response  of  the  rTBC  should  be  due  to  the  positive  response  of  the  electron  density  above  the  F- 
region  maximum. 


3.2  Basic  Approach  of  the  Model  Construction 

The  basic  idea  of  each  global  empirical  TBC  model  which  describes  the  response  to  the 
geomagnetic  activity  is  to  define  a  set  of  analytical  expressions  which  describe  the  most  probable 
TBC  values  at  given  geomagnetic  activity  index,  day  of  the  year,  geographic  location  and  LT.  The 
investigations  on  the  fFl  response  to  the  geomagnetic  activity  (Muhtarov  et  al.,  2002)  indicated 
that  this  is  a  delayed  response  which  can  be  satisfactorily  modeled  by  assuming  that  the 
geomagnetic  influence  is  imposed  on  the  inertial  system  described  by  an  inhomogeneous 
differential  equation  from  a  first  order  (Kutiev  and  Muhtarov,  2003).  In  the  present  study  the  cross¬ 
correlation  analysis  however  indicated  that  the  rTBC  response  to  the  geomagnetic  activity  has  to  be 
presented  by  a  sum  of  two  responses  with  different  time  delay  constants  and  with  different  sign  of 
the  cross-correlation  function.  It  is  known  also  that  during  the  recovery  phase  of  the  ionospheric 
storms  with  geomagnetic  origin  the  ionospheric  reaction  continues  some  time  after  the  geomagnetic 
storm  attenuation.  This  phenomenon  aggravates  the  relationship  between  the  Kp-index  and  the 
ionospheric  anomalies.  In  order  to  resolve  this  problem  Muhtarov  et  al.  (2002)  suggested  an 
approach  for  defining  new  modified  function  of  Kp-index,  based  on  the  time  delay  constant  from  the 
cross-correlation  analysis,  with  its  variations  closely  resembling  those  of  the  relative  foF2.  Having 
in  mind  the  above  mentioned  ideas  Andonov  et  al.  (2011)  constructed  regional  rTBC  model  over 
North  America  and  the  similar  approach  will  be  used  in  the  present  study  as  well. 

If  we  assume  that  the  impact  of  the  geomagnetic  activity  on  the  rTBC  is  accomplished  by  two 
mechanisms  with  different  time  delay  constants  then  the  variability  of  rTBC  can  be  described  as 
follows: 

rrec(r)=  (/,,(A:p,,(r))+/„(Aj)„(()))/,,(ir)/,„{ion)  (12) 

where  the  functions ///LT)  and//o„  (Lon)  represent  the  dependence  of  the  rTBC  response  on  the  LT 
and  longitude  at  equal  other  conditions.  The  parameters  Kpjs  and  Kpn  are  the  modified  with  the  time 
delay  constants  respectively  Ts  and  Ti  values  of  the  Kp-index.  These  modified  parameters  are 
solutions  of  the  equations  shown  below  and  are  obtained  easily  by  a  numerical  integration: 

T,^^^^f^  +  KppXt)=Kp{t)  (13) 

at 

T^‘‘^Prii‘)  +  KpJt)=Kp{t)  (14) 

dt 

The  unknown  functions /t^  and/n  from  (1)  can  be  expressed  by  their  Taylor  series  expansions  while 
the  dependences  on  the  LT  and  longitude  can  be  presented  by  a  Fourier  series  as  follows: 
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Then  the  rTEC  can  be  described  as: 
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We  note  that  the  Fourier  time  series  (third  relation  in  (15))  includes  the  contribution  of  four 
harmonics,  24,  12,  8  and  6  hours,  while  the  Fourier  longitude  series  (fourth  relation  in  (15)) 
includes  the  contribution  of  6  harmonics,  i.e.  the  contribution  of  zonal  waves  with  zonal 
wavenumbers  up  to  6.  It  is  worth  noting  that  the  numbers  of  the  included  components  in  the  Fourier 
expansion  series  are  defined  experimentally.  We  accepted  only  the  contribution  of  the  above 
mentioned  diurnal  components  and  zonal  waves  because  the  addition  of  more  components  does  not 
improve  the  model  error. 

The  next  important  step  is  to  find  a  functional  dependence  between  the  Kp  and  rTEC  in  order  to 
clarify  the  number  of  terms  in  the  Taylor  series  (first  two  relations  in  (15)).  The  most  appropriate 
type  of  the  functional  dependence  is  defined  empirically  by  following  the  approach  described  in 
Andonov  et  al.  (2011).  Some  examples  of  the  empirical  dependences  between  the  Kp  and  rTEC 
calculated  for  different  months  and  geographical  points  which  are  denoted  at  the  plots,  are 
presented  in  Figure  20.  It  is  evident  that  in  all  cases  the  functional  dependence  between  Kp  and 
rTEC  is  close  to  the  cubic  function.  Due  to  this  result  in  the  Taylor  series  only  the  contribution  of 
first  four  terms  are  included.  Then  the  most  probable  values  of  the  coefficients:  C(is,  ocu,  jiu  yu  Tg  and 
Ti  from  (15)  have  to  be  determined.  This  is  a  nonlinear  optimizing  task  that  can  be  solved  by 
applying  the  “trial-and-error”  method  in  a  way  that  the  best  approximation  in  a  sense  of  minimum 
least  squares  deviation  to  be  assured.  In  order  to  solve  the  problem  the  following  steps  are  made:  (i) 
it  is  given  a  range  of  Tg  changes  from  0  to  10  hours  with  a  time  resolution  of  1  hour  and  a  range  of 
Ti  changes  from  11  to  72  hours  with  a  time  resolution  also  of  1  hour;  (ii)  for  each  point  of  the  built 
in  this  way  grid  the  coefficients  Ois,  ocu  j5i  and  y,  are  found  by  using  the  method  of  least  squares  best 
fit,  and  (iii)  the  coefficients  otis,  On,  fit,  ju  Tg  and  Ti  at  which  the  best  approximation  (in  a  sense  of 
minimum  least  squares  deviation)  is  obtained  are  accepted  as  optimal  coefficients  for  the  model 
rTEC  described  by  (16). 

In  the  present  study  we  accepted:  (i)  longitude  and  UT  as  independent  variable  quantities;  the 
conversion  to  LT  is  a  simple  procedure,  and  (ii)  at  each  modip  latitude  a  separate  model,  described 
by  (5),  is  constructed;  the  values  of  the  model  rTEC  which  do  not  belong  to  the  5°  modip  grid  are 
obtained  by  an  interpolation  procedure  that  will  be  described  later.  The  latter  is  done  because  if  a 
latitudinal  approximation  is  used  first  the  number  of  model  constants  will  increase  and  second  an 
additional  error  will  be  introduced.  The  rTEC  model  described  by  (16)  contains  820  constants  and 
they  are  determined  by  least  squares  fitting  techniques. 
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Figure  20  The  empirical  dependence  between  the  Kp  and  rTEC  calculated  for  different  months  and  geographical  points 
noted  at  the  plots. 


3.3  rTEC  Model  Results 


In  order  to  demonstrate  how  the  model  is  able  to  describe  the  rTEC  response  to  geomagnetic 
activity  two  geomagnetic  storms  observed  at  different  seasons  and  solar  activity  conditions  will  be 
presented. 


Figure  21a  Temporal  variability  of  the  Kp-index  during  the  geomagnetic  storm  in  April  5-13,  2000. 
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Figure  6b  Comparison  between  the  model  (right  column  of  plots)  and  observed  (left  column  of  plots)  rTEC  longitude- 
hour  cross-sections  for  the  considered  geomagnetic  storm  at  different  modip  latitudes  noted  above  the  plots. 
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Figure  21a  shows  the  temporal  variability  of  the  Kp-index  during  the  geomagnetic  storm  in  April  5- 
13,  2000,  at  high  solar  activity.  The  rapid  increase  of  the  Kp-index  starts  at  around  12  UT  on  6 
April  and  reaches  the  largest  values  (close  to  9)  at  midnight  and  early  hours  on  7  April;  then  it 
decreases  to  the  undisturbed  levels.  Figure  21b  presents  the  comparison  between  the  model  (right 
column  of  plots)  and  observed  (left  column  of  plots)  rTBC  longitude-hour  cross-sections  for  the 
considered  period  of  time,  April  5-13,  at  different  modip  latitudes  which  are  noted  above  the  plots. 
We  clarify  that  the  longitude  is  presented  by  numbers  define  from  longitude/15°  while  the  time  is  in 
hours  and  it  starts  on  01  April  2000  at  00  UT.  In  order  to  facilitate  the  comparison  between  the 
model  and  data  results  the  color  scales  are  the  same  at  the  same  modip  latitudes.  However  as  the 
rTBC  response  strongly  depends  on  the  modip  latitude  the  color  scales  are  different  at  different 
modip  latitudes.  The  careful  comparison  between  model  and  data  plots  reveals  that  the  model 
overall  well  reproduces  the  real  situation.  It  however  underestimates  the  positive  rTBC  response  at 
high  latitudes  in  the  second  half  of  April  6  at  70°N  i.e.  almost  simultaneously  with  the  largest  values 
of  Kp-index  and  one  day  later  positive  response  at  70°S.  The  model  quite  well  reproduces  the 
negative  rTBC  response  in  the  second  half  of  April  6  at  70°S  when  Kp-index  rapidly  increases.  At 
middle  latitudes,  +40°,  the  model  positive  rTEC  response  is  also  slightly  underestimated  in  the 
second  half  of  April  6.  The  model  rTBC  response  at  tropical  and  equatorial  latitudes  comparatively 
well  reproduces  the  data.  It  is  worth  noting  that  at  all  modip  latitudes  the  model  very  well  describes 
the  longitude  variabilities  of  the  rTBC  response.  Most  probably  this  is  due  to  the  large  number  of 
the  zonal  waves  included  in  the  model. 


Figure  22a  Temporal  variability  of  the  Kp-index  during  the  geomagnetic  storm  in  November  6-13,  2004. 

Figure  22a  presents  the  temporal  variability  of  the  Kp-index  during  the  geomagnetic  storm  in 
November  6-13,  2004,  i.e.  during  the  declining  phase  of  the  solar  cycle.  This  is  more  complex 
geomagnetic  storm  with  two  Kp-index  amplifications  which  are  far  from  each  other  of  about  two 
days.  Figure  22b  presents  the  comparison  between  the  model  (right  column  of  plots)  and  observed 
(left  column  of  plots)  rTBC  longitude-hour  cross-sections  for  the  period  of  time,  November  6-13, 
again  at  different  modip  latitudes  arranged  in  the  same  way  as  those  in  Figure  21b.  In  this  case  the 
rTBC  responses  to  both  Kp-index  amplifications  have  to  be  considered.  The  model  comparatively 
well  reproduces  the  rTBC  response  related  to  the  first  Kp-index  amplification  at  high  latitudes,  +70°, 
but  underestimate  that  connected  with  the  second  Kp-index  amplification.  The  temporal  and 
longitudinal  variability  of  the  rTBC  response  at  middle  latitudes,  +40°,  is  very  well  recreated;  only 
the  positive  rTBC  response  on  November  8  is  slightly  underestimated.  Some  slight  underestimation 
of  the  rTBC  response  can  be  also  distinguished  at  tropical  and  equatorial  latitudes. 

The  above  shown  examples  of  two  geomagnetic  storms  clearly  indicate  that  the  global  empirical 
rTBC  model  describes  quite  well  the  ionospheric  response  to  the  geomagnetic  activity  at  different 
solar  cycle  conditions.  Each  empirical  model  needs  to  be  statistically  evaluated.  The  main  statistics 
based  on  the  entire  data  set  are  presented  here  through  the  mean  (systematic)  error  (MB),  root  mean 
squares  error  (RMSE)  and  the  standard  deviation  error  (STDE)  usually  accepted  as  the  basic  error 
characteristics  of  each  model.  In  order  to  evaluate  the  main  statistics  of  the  model  first  the  model 
TEC,  where  TECmod  =  TECmed(l+rTEC),  i.e.  this  is  a  corrected  15-day  median  TEC  with  the  model 
rTBC,  is  calculated.  The  following  errors  for  the  considered  period  of  time  (1  January  1999-31 
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Figure  22b  Comparison  between  the  model  (right  column  of  plots)  and  observed  (left  column  of  plots)  rTEC  longitude- 
hour  cross-sections  for  the  considered  geomagnetic  storm  at  different  modip  latitudes  noted  above  the  plots. 
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December  2011)  are  obtained:  M£'=-0.204  TBCU,  i.e.  the  model  fits  to  the  CODB  TBC  data  with 
small  negative  bias;  then  the  RMSE  and  STDE  have  very  close  values,  i.e.  RMSE=4.592  and 
STDE=4.588. 

The  overall  statistics  of  the  model  error  can  be  defined  more  precisely  by  showing  the  dependence 
of  RMSE  on  modip  latitude  and  months.  Figure  23  shows  a  modip  latitude-month  cross-section  of 
the  model  RMSE  calculated  for  the  entire  (January  1999  -  December  2011)  period  of  time.  The 
largest  RMSE,  reaching  7.5  TBCU,  are  observed  at  low  latitudes  where  the  equatorial  ionospheric 
anomaly  is  developed;  the  crests  are  situated  at  around  +30°  modip  latitude.  Some  amplification  of 
the  RMSE  can  be  noticed  also  around  70°S  during  the  equinoxes;  it  reaches  around  4  TBCU.  The 
largest  errors  are  found  mainly  in  the  equinoxes  and  this  could  be  due  to  both  semiannual  variability 
of  the  ionosphere  and  semiannual  variability  of  the  geomagnetic  activity.  We  calculated  also  the 
RMSE  for  each  month  of  the  entire  time  interval  (not  shown  result);  as  it  is  expected  in  advance  the 
RMSE  has  larger  values  during  high  solar  activity. 
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Figure  23  Modip  latitude-month  cross-section  of  the  model  RMSE  calculated  for  the  entire  (January  1999  -  December 
2011)  time  interval. 


3.4  Global  TEC  Maps 

The  basic  aim  of  each  global  TBC  model  used  for  short-term  prediction  is  to  construct  the  global 
distribution  of  the  TBC,  i.e.  to  obtain  global  TBC  maps,  at  given  geomagnetic  activity,  day  of  the 
year  and  UT.  The  presented  here  rTEC  model  predicts  the  correction  to  the  15-day  median  values 
for  each  hour  of  the  prediction  period.  The  TEC  value  at  a  given  hour  is  actually  a  corrected  15-day 
median  TEC  with  the  model  rTEC.  Then  the  TEC  maps  are  constructed  by  interpolation  of  the  TEC 
values  from  the  used  grid  with  a  5°  step  in  modip  latitude  and  15°  in  longitude.  The  interpolation 
between  obtained  TEC  values  is  done  again  by  using  Inverse  Distance  Method  (Shepard,  D.,  1968). 
Then  the  modip  frame  is  converted  to  geographical  one.  The  TEC  values  assigned  to  both  poles  are 
found  by  interpolation  between  the  known  from  the  model  points  which  have  the  highest  northern 
and  southern  latitudes.  The  model  maps  are  arrayed  in  terms  of  the  coordinate  system  of 
geographical  latitude  from  -90°  to  90°  at  each  5°  and  longitude  from  -180°  to  180°  at  each  5°. 

Figure  24  presents  a  comparison  between  the  model  (right  column  of  plot)  and  CODE  TEC  maps 
for  November  8,  2004  geomagnetic  storm  at  OOUT  (upper  row  of  plots),  06UT  (second  row  of 
plots),  12UT  (third  row  of  plots)  and  18UT  (bottom  row  of  plots).  In  order  to  facilitate  the 
comparison  the  color  scales  at  each  pair  of  TEC  maps  (at  one  and  the  same  UT)  are  the  same.  Due 
to  the  dependence  of  the  TEC  on  the  UT  the  color  scales  are  not  the  same  at  different  UT.  In 
general  there  is  significant  similarity  between  the  CODE  and  model  TEC  maps  but  there  are 
differences  as  well.  At  00  UT  and  06  UT  for  example,  the  equatorial  anomaly  is  underestimated  by 
the  model  but  the  hemispheric  asymmetry  is  quite  well  reproduced.  At  12  UT  and  18  UT  the  model 
densest  part  of  the  ionosphere  is  closer  to  the  observations;  some  longitude  structures,  as  for 
example  the  Weddell  Sea  Anomaly  in  the  SH  at  18  UT,  can  be  reproduced  well. 
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Figure  24  Comparison  between  the  model  (right  column  of  plot)  and  CODE  TEC  maps  for  November  8,  2004 
geomagnetic  storm  at  OOUT  (upper  row  of  plots),  06UT  (second  row  of  plots),  12UT  (third  row  of  plots)  and  18UT 
(bottom  row  of  plots).  The  modip  latitude  is  also  marked  by  white  line. 

The  above  presented  comparison  between  the  model  and  CODE  TEC  maps  revealed  good  similarity. 
This  means  that  the  constructed  global  rTEC  model  response  to  the  geomagnetic  activity  can  be 
successfully  used  for  short-term  (24  hours  ahead)  prediction. 

3.5  Concluding  comments  and  on-line  implementation  of  the  global  model  of  TEC 
response  to  geomagnetic  activity 

Section  3  presented  a  global  empirical  TEC  model  response  to  geomagnetic  activity  described  by 
the  Kp-index  (Mukhtarov  et  al.,  2013c).  It  describes  the  geomagnetically  forced  changes  of  the  TEC 
assuming  that  these  changes  at  a  given  modip  latitude  depend  mainly  on  Kp-index,  LT  and 
longitude.  The  geomagnetic  changes  are  expressed  by  the  relative  deviation  of  TEC  from  its  15-day 
median,  noted  as  rTEC.  Therefore  this  model  predicts  the  correction  to  the  15-day  median  values, 
rTEC,  for  each  hour  of  the  prediction  period.  The  model  offers  TEC  maps  which  depend  on 
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geographic  coordinates  (5°x5°  in  latitude  and  longitude)  and  UT  at  given  geomagnetic  activity  and 
day  of  the  year. 

The  rTBC  model  contains  820  constants  and  they  are  determined  by  least  squares  fitting  techniques. 
The  comparison  between  the  model  and  CODE  TEC  data  at  different  solar  cycle  and  season 
conditions  demonstrated  high  degree  of  similarity.  The  very  small  systematic  error  (-0.204)  and  low 
RMSE  (4.592  TECU)  of  the  model  characterized  it  as  useful  tool  for  describing  the  ionospheric 
TEC  response  to  geomagnetic  storms. 

This  model  can  be  used  for  short-term  prediction.  For  this  purpose  at  a  given  day  of  the  year, 
geographic  location  and  UT  the  model  needs  as  input  parameter  only  the  predicted  Kp-index.  This 
is  a  possible  task  because  there  are  available  models  which  predict  the  geomagnetic  activity  with  a 
reliable  accuracy.  An  example  of  such  model  is  a  MAK  model  described  by  Andonov  et  al.  (2004). 
It  provides  online  prediction  of  the  .^T^-index  and  is  implemented  on  the  web  site: 

http://www.geophys.bas.bg/kp_for/kp_mod_bg-php- 
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Figure  25  The  block  chart  of  the  global  TEC  model  for  short-time  prediction 
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The  online  software  of  the  model  has  been  already  implemented  on  the  website 
(http://www.geophys.bas.bg/tec/page/tec_hourly_en.htm)  and  it  is  also  still  in  testing  mode.  Figure 
25  shows  its  block-chart.  It  will  be  officially  implemented  at  the  institute  website  when  the  paper 
Mukhtarov  et  al.  (2013c)  is  published.  The  short-term  TBC  prediction  particularly  at  strong 
geomagnetic  storms  will  improve  significantly  the  accuracy  of  the  geodetic  and  navigation  data 
which  have  increasing  importance  in  resolving  both  scientific  and  practical  tasks. 


4.  On-line  demonstration  of  the  TEC  model  results 

The  two  TBC  models  have  been  officially  implemented  at  the  institute  website  since  October  2013 
when  all  three  papers  devoted  to  them  were  published  in  J.  Geophys.  Res.  -  Space  Physics.  They 
can  be  found  at:  http://www.geophys.bas.bg/tec/page/tec.html.  At  this  web  side  besides  the  short¬ 
term  (24  hours  ahead)  and  long-term  (a  month  ahead)  prediction  TBC  maps  produced  by  the  TBC 
model  response  to  geomagnetic  activity  and  the  background  TBC  model  respectively  there  are  also 
two  example  TBC  maps,  showing  as  a  movie,  demonstrating  the  strong  sides  of  both  models.  The 
hourly  median  TBC  maps  for  January  01,  2012  are  shown  where  the  development  of  the  Weddell 
Sea  Anomaly  is  clearly  reproduced  by  the  background  TBC  model.  The  geomagnetically  forced 
TBC  anomalies,  expressed  by  the  relative  deviation  of  TBC  from  its  15-day  median  (rTEC)  for 
October  29  -  November  01,  2003  demonstrate  the  global  TBC  response  to  the  famous  Halloween 
geomagnetic  storm.  A  detailed  comparison  between  the  data  and  model  can  be  seen  by  a  link  there. 

It  has  been  already  mentioned  that  the  TBC  model  for  short-term  prediction  needs  as  input 
parameter  only  the  predicted  Kp-index.  In  our  case  it  is  taken  from  the  MAK  model  (Andonov  et  al., 
2004),  which  provides  online  prediction  of  the  ^^-index  for  the  next  six  hours  from  the  current  time. 
Figure  26  shows  how  the  MAK  model  works  at  the  web  site  and  how  one  of  the  last  geomagnetic 
storm,  19  February  2014,  is  described  there. 
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Figure  26  Geomagnetic  Kp-index  forecast  based  on  the  MAK  model. 

The  strongest  geomagnetic  storm  during  the  last  20  years  is  the  Halloween  one,  29-31  October  2003. 
The  website  presents  both  hourly  maps  of  the  models  and  those  of  the  CODB  data  for  the  period  of 
time  29  October-01  November  2003  when  the  geomagnetic  storm  takes  place.  Two  types  of  maps 
are  shown  on  the  web  site:  rTEC  and  TEC.  Figure  27  presents  the  comparison  between  the  rTEC 
and  TEC  maps  on  30  October  2003  at  22  UT  (when  the  second  Kp-index  peak  is  seen)  produced  by 
the  CODE  data  (left  plot)  and  model  results  (right  plot). 
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Figure  27  (upper  plot)  The  Kp-index  for  29  Oct-01  Nov  2003;  comparison  between  the  rTEC  (upper  row)  and  TEC 
maps  (bottom  row)  between  the  CODE  data  (left  column)  and  model  results  (right  column). 


The  careful  inspection  of  the  two  couples  of  maps  shows  a  high  degree  of  similarity;  the  model 
reacts  to  both  positive  and  negative  TBC  anomalies.  This  indicates  that  the  TBC  model  response  to 
geomagnetic  activity  captures  very  well  all:  FT,  modip  latitude/longitude  and  seasonal  dependences 
of  the  TBC  response  to  strong  geomagnetic  storms  like  the  Halloween  one.  The  quantitative 
comparison  however  shows  some  underestimation  of  the  magnitude  of  the  TBC  anomalies.  Of 
course,  this  is  a  consequence  of  the  fact  that  all  empirical  models  predict  an  average  response  based 
on  the  data  used  for  constructing  the  model.  We  have  carefully  assessed  the  ability  of  the  model  to 
reproduce  and  predict  the  TBC  response  at  different  conditions.  It  has  been  found  that  the  model  is 
able  to  capture  well  the  TBC  response  even  to  moderate  and  minor  geomagnetic  storms  when 
however  the  solar  radio  flux  variability  is  low,  or  particularly  when  the  ~27-day  (~13.5-day  also) 
oscillations  are  absent.  This  is  related  to  the  fact  that  the  model  predicts  the  correction  to  the  15-day 
median  values,  but  these  15-day  medians  significantly  suppress  such  oscillations. 

Figure  28  presents  an  example  of  how  the  web  site  for  both  TBC  models  looks  like  for  26  March 
2014.  After  the  representation  of  the  Weddell  Sea  Anomaly  (upper  left  plot)  and  rTBC  maps  for 
Halloween  geomagnetic  storm  (upper  right  plot),  which  demonstrate  the  ability  of  both  models  to 
reproduce  different  ionospheric  events,  then  the  short-term  (for  27  March,  middle  plot)  and  long¬ 
term  (for  26  April,  bottom  plot)  prediction  TBC  maps  are  shown. 
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National  Institute  of  Geophysics,  Geodesy  and  Geography 
Global  empirical  TEC  models  based  on  the  CODE  TEC  data 

The  ionospheric  total  electron  content  (TEC)  is  a  key  parameter  for  describing  the  impact  of  the  ionosphere  on  the  propagation  of  radio  signals. 
Understanding  the  variability  of  TEC  is  crucial  for  the  operation  of  many  applications,  as  navigation  satellite  systems  like  Global  Positioning  System 
(GPS),  Global  Navigation  Satellite  System  (GLONASS),  and  the  future  Galileo  system. 

Two  global  empirical  TEC  models  are  built  by  using  the  Center  for  Orbit  Determination  of  Europe  (CODE)  TEC  data  (ftp://ftp.unibe.ch/alub/CODE/) 
for  full  13  years,  1999-2011.  The  background  TEC  model  (for  details  see  the  two  PDF  files,  [PDFl],  [PDF2])  offers  31-day  running  median  TEC  maps 
for  a  month  ahead  at  a  given  UT.  The  TEC  model  response  to  geomagnetic  activity  (for  details  see  the  PDF  file)  presents  a  current  and  24  hours 
ahead  TEC  maps. 

The  work  on  the  TEC  models  was  supported  by  the  European  Office  of  Aerospace  Research  and  Development  (EOARD),  Air  Force  Office  of  Scientific 
Research,  Air  Force  Material  Command,  USAF,  under  grant  number  FA8655-12-1-2057  to  D.  Pancheva  (for  contact). 

1.  Global  empirical  background  TEC  model  2.  Global  empirical  model  of  TEC  response  to  geomagnetic 

activity 


The  hourly  median  TEC  maps  for  January  01,  2012  are  shown  where  The  geomagnetically  forced  TEC  anomalies,  expressed  by  the  relative 
the  development  of  the  Weddell  Sea  Anomaly  can  be  clearly  deviation  of  TEC  from  its  15-day  median  (rTEC)  for  October  29  -  November 
reproduced  by  the  background  TEC  model.  01,  2003  demonstrate  the  global  TEC  response  to  the  famous  Halloween 

geomagnetic  storm.  A  detailed  comparison  between  the  data  and  model 
see  HERE 
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Figure  28  Global  empirical  TEC  models  based  on  the  CODE  TEC  data 
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Very  recently  a  mobile  version  of  both  TBC  models  have  been  implemented  and  now  the  TBC 
model  predictions  have  be  seen  by  using  tablets  or  smart  phones.  This  was  done  in  order  to  support 
advertizing  the  models. 


5.  An  attempt  to  mitigate  the  global  background  TEC  model’s  error  by  using 
regularly  arriving  new  TEC  data  -  a  hybrid  TEC  model 

The  attempt  for  improving  the  quality  of  the  long-term  TBC  forecast  is  based  on  the  regularly 
arriving  fresh  CODB  TBC  data  and  the  autocorrelation  prediction  of  the  error  and  the  respective 
correction  of  the  background  model  with  the  predicted  error.  Because  this  approach  needs  fresh  data, 
i.e.  new  data  which  have  not  been  used  for  constructing  the  model,  it  is  called  a  hybrid  method. 
The  CODB  TBC  data,  used  for  constructing  the  background  TBC  model,  arrive  daily  with  a  delay  of 
a  few  days  toward  the  current  date.  As  this  model  predicts  the  31 -day  running  median  TEC  values 
then  reliable  (measured)  TEC  data  for  a  past  period  are  available  only  for  about  15  days  from  the 
date  at  which  the  prediction  is  made.  If  the  time  distance  between  the  dates  for  which  the  prediction 
is  made  to  that  for  which  the  reliable  data  for  calculating  the  31 -day  running  median  TEC  value  are 
available  is  denoted  as  ojfset,  then  the  offset  will  be  equal  or  larger  than  15  days  as  it  is  composed 
by  number  of  days  between  the  current  date  and  the  date  for  which  the  prediction  is  made  plus  15 
days  before  the  current  date. 

The  suggested  method  will  be  demonstrated  on  the  archival  data.  The  main  idea  is  based  on  the 
assumption  that  the  time  series  of  the  model  residuals  are  not  composed  only  of  random  component 
but  have  also  some  time  dependent  component.  Figure  29  shows  the  temporal  variability  of  the 
model  residuals  calculated  for  points  with  longitude  of  00°E  and  modip  latitudes  of  40°N  (upper 
plot),  00°  (middle  plot),  and  40°S  (bottom  plot)  for  00  LT.  The  careful  inspection  of  the  plots 
reveals  that  the  model  residuals  are  correlated,  i.e.  each  value  is  related  to  the  neighbour  ones.  This 
circumstance  provides  an  opportunity  for  predicting  the  model  residuals  in  future  moments  on  the 
base  of  the  actual  residuals  in  the  past  which  are  determined  by  comparison  of  the  model  values 
with  the  observations.  The  determination  of  the  predicted  residuals  could  provide  a  possibility  for 
correcting  the  model  values  for  some  future  moments. 


Figure  29  Temporal  variability  of  the  model  residuals  calculated  for  points  with  longitude  of  00°  and  modip  latitudes  of 

00°  (upper  plot),  40°N  (left  plot)  and  40°S  (right  plot)  for  00  LT. 
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The  autocorrelation  method  for  forecasting  the  error  is  based  on  the  well-known  theorem  of 
Wiener-Hopf.  For  a  stationary  random  process  x(t),  for  which  the  mathematical  expectation, 
variance  and  autocorrelation  function  are  known  quantities  then  an  unknown  value  of  the  process 
x(t)  at  future  moment  tf  can  be  represented  at  minimum  squares  deviation  as  a  linear  combination  of 
some  sample  of  known  values  of  the  process  x(t)  at  the  moments:  ti,  t2, ...,  4.  Then: 

x(tf)  =  x  +  ^Pi{x{ti)-x)  (17) 


The  coefficients  Pi  from  (3)  are  a  solution  of  the  system  of  equations: 


=  k  =  \,2,...,n  (18) 

p=i 


The  normalized  (by  dispersion)  autocorrelation  function  of  the  process  x{t),  denoted  by  p,  is  present 
in  the  system  of  equations  (18)  and  it  is  a  function  of  time  lag.  The  right  part  of  the  system  of 
equations  (18)  contains  the  autocorrelations  at  time  lags  between  the  moment  for  prediction  and  the 
moments  where  there  are  known  values.  Hence  the  deviation  of  the  prediction  from  the 
observations,  based  on  this  method,  depends  on  the  values  of  the  autocorrelation  coefficients  at  time 
lags  corresponding  to  the  distances  between  the  moment  for  prediction  and  the  moments  for  which 
the  values  of  the  process  x{t)  ate  known. 


The  simulations  of  the  autocorrelation  correction  of  the  model’s  error  will  be  demonstrated  on  the 
TEC  data  for  the  time  period  2007-2012.  The  autocorrelation  function  of  the  deviations  (residuals 
or  errors)  is  calculated  by  using  the  model  residuals  for  the  period  of  time  1999-2006. 
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Figure  30  Normalized  autocorrelation  functions  calculated  for  equator  (uppermost  plot)  and  at  modip  latitudes  of:  40°N 
(left  upper  plot);  40°S  (right  upper  plot);  80°N  (left  bottom  plot)  and  80°S  (right  bottom  plot). 
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Figure  30  shows  the  normalized  autocorrelation  functions  calculated  for  equator  (uppermost  plot) 
and  at  modip  latitudes  of:  40°N  (left  upper  plot);  40°S  (right  upper  plot);  80°N  (left  bottom  plot)  and 
80°S  (right  bottom  plot).  The  inspection  of  the  autocorrelation  functions  reveals  that  a  higher 
correlation  of  the  model  residuals  is  observed  at  lower  modip  latitudes  than  those  at  high  modip 
latitudes.  Values  of  the  correlation  coefficients  around  and  above  0.5  are  seen  at  time  lags  up  to  60 
days,  i.e.  an  effective  correction  of  the  error  is  possible  to  occur  if  the  offset  is  not  more  than  two 
months. 

The  next  step  is  to  check  how  the  length  of  the  offset  affects  the  correction  approach.  The  validation 
procedure  is  done  at  offsets  of  15,  30,  45,  60,  75  and  90  days.  It  actually  presents  a  simulation  of  the 
predicted  error  in  real  conditions.  The  correction  of  the  global  background  TBC  values  by  the 
autocorrelation  method  for  error  prediction  is  demonstrated  in  Figure  31.  It  shows  a  comparison 
between  the  observed  TBC  (thin  solid  line),  not  corrected  model  TBC  (thick  solid  line)  and 
corrected  model  TBC  (dashed  line)  at  offsets  of:  15  days  (left  upper  plot),  30  days  (right  upper  plot), 
45  days  (left  bottom  plot)  and  60  days  (right  bottom  plot)  calculated  for  30°N  and  0°B  at  12  UT.  It 
is  seen  that  if  the  offset  is  15  days  (i.e.  to  predict  the  error  at  the  current  day)  the  corrected  model 


Figure  31  Comparison  between  the  observed  TEC  data  (thin  solid  line),  not  corrected  model  TEC  (thick  solid  line)  and 
corrected  model  TEC  (dashed  line)  at  offsets  of:  15  days  (left  upper  plot),  30  days  (right  upper  plot),  45  days  (left 
bottom  plot)  and  60  days  (right  bottom  plot)  calculated  for  30°N  and  0°E  at  12  UT. 
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TBC  almost  coincides  with  the  TBC  data.  With  increasing  the  offset  the  difference  between  the 
corrected  model  TBC  and  the  TBC  data  increases  as  the  corrected  model  TBC  values  approach  to 
values  of  the  not  corrected  ones.  The  presented  in  Figure  31  examples  clearly  demonstrate  the 
efficiency  of  the  correction  depending  on  the  given  offset. 

The  advantage  of  the  hybrid  model  for  long-term  prediction  with  respect  to  the  built  by  Muklitarov 
et  al.  (2013a)  background  TEC  model  is  illustrated  by  comparing  of  their  RMSE.  Figure  32  shows 
the  comparison  between  the  RMSE  of  the  background  TEC  model  (thick  solid  line)  and  that  of  the 
corrected  TEC  model  (thin  solid  line  with  symbols)  calculated  for  different  offsets  for  the  period  of 
time  2007-2012.  The  found  dependence  of  the  RMSE  on  the  offset  support  the  earlier  suggestion 
that  the  correction  is  really  effective  if  the  error  prediction  is  made  for  a  date  with  a  distance  up  to 
60  days  from  the  date  with  real  data.  If  we  consider  real  conditions  a  prediction  of  a  month  ahead 
means  that  the  offset  is  45  days.  In  this  case,  the  RMSE  decreases  from  3.2  TECU  to  2.76  TECU.  If 
we  make  a  prediction  for  the  current  month  then  the  offset  is  15  days  and  the  RMSE  falls  down  to 
1.52  TECU. 


15  30  45  60  75  90 

Offset  (days) 


Figure  32  Comparison  between  the  RMSE  of  the  not  corrected  background  TEC  model  (thick  solid  line)  and  corrected 
global  TEC  model  (thin  solid  line  with  symbols)  calculated  for  different  offsets  for  the  period  for  time  2007-2012. 

The  procedure  for  the  correction  of  the  model  residuals  is  fully  applicable  to  the  real  prediction  of 
the  global  median  TEC.  When  it  is  necessary  to  make  TEC  prediction  for  a  date  with  large  offset 
then  the  effectiveness  of  the  correction  asymptotically  approaches  to  zero.  This  is  due  to  the  fact 
that  that  the  right  side  of  the  system  of  equations  (18)  is  inclined  to  zero  and  then  the  coefficients  Pi 
from  (17)  become  zero.  If  this  happens  however  the  quality  of  the  median  global  TEC  prediction 
does  not  get  worse.  The  presented  in  this  part  a  hybrid  method  is  described  in  detail  by  Mukhtarov 
et  al.  (2014). 

In  conclusion  it  is  worth  noting  that  the  hybrid  model  is  fully  applicable  to  the  real  prediction  of  the 
global  median  TEC.  The  hybrid  model  however  is  not  an  autonomous  ionospheric  TEC  model  and 
due  to  this  it  is  not  useful  for  numerous  single  frequency  GPS  applications  that  need  additional 
information  (as  an  autonomous  ionospheric  TEC  model)  for  mitigating  the  ionospheric  propagation 
error. 


6.  An  attempt  to  assess  the  stratospheric  impact  on  the  TEC  variability  in  winter 

The  last  task  proposed  to  be  investigated  in  this  project  is  related  to  the  ionospheric  response  to 
anomalous  atmospheric  phenomena,  like  sudden  stratospheric  warming  (SSW)  events  (Pancheva 
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and  Mukhtarov,  2011b;  Jin  et  al.,  2013)  that  has  been  recently  obtained.  This  finding  has  directed 
the  attention  of  the  researchers  to  the  relationship  between  the  TEC  and  stratospheric 
meteorological  field  (as  temperature,  geopotential  height,  neutral  wind  etc.)  variabilities.  Due  to  this 
the  next  step  is  to  conduct  a  detailed  research  on  this  ionosphere-atmosphere  relationship, 
particularly  strong  in  winter,  and  to  present  an  idea  about  a  TEC  model  where  the  stratospheric 
impact  is  included.  As  the  stratospheric  meteorological  fields  are  predicted  with  a  few  days  ahead  it 
is  worth  making  an  attempt  to  add  also  the  stratospheric  effect  on  the  TEC  variability  in  a  model  for 
short-time  TEC  prediction.  We  expect  that  some  improvement  of  the  TEC  model,  particularly 
during  the  SSW  events,  could  be  possible.  In  this  part  we  will  present  only  some  research  results 
that  later  could  be  incorporated  for  improving  the  global  TEC  model  for  short-term  prediction. 

The  numerical  simulations  performed  by  Liu  and  Roble  (2002)  for  the  first  time  reported  that  the 
stratospheric  warming  is  associated  with  the  lower  thermospheric  warmings.  Later  this  result  was 
supported  by  the  satellite  measurements  of  the  MIPAS/Envisat  which  indicated  that  the  high 
latitude  thermosphere  between  120  and  140  km  is  really  warmer  during  the  major  SSW  in  January 
2009  (Funke  et  al.,  2010).  The  temperature  changes  in  the  lower  thermosphere,  i.e.  in  the  dynamo 
region,  lead  however  to  wind  changes  which  through  dynamo  effect  could  have  impact  on  the 
ionospheric  plasma  redistributions.  This  simple  physical  relationship  was  used  by  Pancheva  and 
Mukhtarov  (2011b)  who  for  the  first  time  reported  a  significant  decrease  of  both  mean  electron 
density  and  the  amplitude  of  the  diurnal  migrating  component  in  the  low  latitude  ionosphere  during 
the  major  SSW  in  January  2009.  To  explain  the  observations  seen  in  the  satellite  FORMOSAT- 
3/COSMIC  electron  density  the  author  suggested  a  mechanism  analogous  to  the  so  called 
“disturbed  dynamo”,  but  generated  by  the  SSW  event  not  by  geomagnetic  storms.  Later  the 
validity  of  this  mechanism  was  supported  by  numerical  simulations  with  the  GSM  TIP  model 
reported  by  Klimenko  et  al.  (2012). 

The  above  results  directed  our  attention  to  the  temperature  as  a  meteorological  parameter  that  has  to 
be  used  for  searching  correlation  with  the  ionospheric  TEC.  For  this  purpose  the  satellite  MLS-Aura 
temperature  measurements  and  the  CODE  TEC  data  have  been  used.  Before  calculating  the  cross¬ 
correlation  functions  the  seasonal  courses  were  removed.  This  is  done  by  using  31 -day  running 
means  which  were  removed  from  the  original  data.  The  preliminary  investigations  concerning  the 
choice  of  the  most  representative  altitude  and  latitude  for  the  stratospheric  temperature  revealed  that 
the  temperature  at  the  altitude  around  40  km  and  latitude  of  60°N  describes  the  most  typical  winter 
conditions.  Due  to  this  the  cross-correlation  with  the  global  TEC  has  been  done  with  this 
temperature  for  the  boreal  winter  (October-March).  The  results  are  shown  in  Figure  33. 
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Figure  33  (left  plot)  Latitude  structure  of  the  cross-correlation  functions  between  the  ZM  TEC  and  ZM  temperature  at 
latitude  of  60°N  and  altitude  of  ~40  km;  the  same  as  the  left  plot  but  the  cross-correlation  is  with  the  DW 1  TEC. 


The  cross-correlation  results  support  the  observations  reported  by  Pancheva  and  Mukhtarov  (2011b) 
and  the  largest  cross-correlation  coefficients  for  ZM  TEC  and  DW  1  TEC  are  respectively  -0.30  and 
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-0.28.  According  to  the  “Fisher’s  Z-transformation”  test  the  obtained  coefficients  are  above  95% 
confidence  level,  i.e.  they  are  significant. 


The  main  effect  on  the  ionosphere  however  has  the  solar  radiation  represented  by  it  proxy  solar 
radio  flux  F10.7.  The  geomagnetic  effect  is  excluded  because  we  have  already  built  a  global  TEC 
model  response  to  geomagnetic  activity  and  because  the  effect  is  not  linear.  Then  the  basic  aim  of 
this  study  is  to  find  an  approximate  quantitative  relationship  between  both  effects  on  the 
ionospheric  TEC,  i.e.  that  of  stratospheric  temperature  and  F10.7,  during  period  of  low  and  middle 
solar  activity.  For  this  purpose  we  investigate  the  period  of  time  2005-2010.  Figure  34  shows  F10.7 
for  the  period  of  time  2005-2010  (upper  plot)  and  the  latitude  structures  of  the  cross-correlation 
functions  of  F10.7  with  ZM  TEC  (left  plot)  and  DWl  TEC  (right  plot).  The  cross-correlations 
indicate  the  known  positive  relation  between  FI 0.7  and  TEC  and  also  the  impact  of  27-day 
oscillations. 
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Figure  34  (a)  Daily  values  of  F10.7  for  the  period  of  time  2005-2010;  (b)  Latitude  structure  of  the  cross-correlation 
functions  of  F10.7  with  ZM  TEC  (left  plot)  and  DWl  TEC  (right  plot). 


The  above  obtained  cross-correlation  results  have  presented  evidences  for  creating  a  linear 
regression  TEC  model  where  the  impact  of  the  stratospheric  temperature  at  altitude  of  ~40  km  and 
latitude  of  60°N  and  the  day-to-day  variability  of  FI 0.7  on  the  main  TEC  decomposition 
components,  i.e.  ZM  and  DWl,  is  included.  The  model  can  be  described  by: 

TEC ^ {mlcit,dciy)  —  {tnlcit)T^{day  )-l-  CC2MS01  {ftil£it)FlOl{day  —  )"^ ^jmoo 

TEC Dw  ipilcit,  day )  —  {pilatyTp^  {day  —  Tj^^ )  -I-  {nilat^F  101  {day  —  Tp^i^p )  -I- 

The  regression  coefficients  and  the  respective  time  lags  are  defined  by  least  squares  best  fit 
approach.  Figure  35  (a)  presents  modip  latitude  dependence  of  the  regression  coefficients  which 
describe  ZM  TEC  (left  plot)  and  DW  1  TEC  (right  plot),  while  (b)  shows  modip  latitude  dependence 
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of  the  model  RMSB  (in  TBCU)  calculated  for  ZM  TBC  (purple)  and  DWl  RFC  (magenta).  The 
results  indicated  that  this  model  is  valid  for  the  tropical  region  (+30  modip)  where  the  equatorial 
anomaly  is  located  and  where  the  SSW  response  found  by  Pancheva  and  Mukhtarov  (2011b)  was 
the  largest. 

(a) 


(b) 


Figure  35  (a)  Modip  latitude  dependence  of  the  regression  coefficients  which  describe  ZM  TEC  (left  plot)  and  DW 1 
TEC  (right  plot);  (b)  Modip  latitude  dependence  of  the  model  RMSE  (in  TECU)  calculated  for  ZM  TEC  (purple)  and 

DWl  REC  (magenta). 
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Figure  36  Modip  latitude  dependence  of  the  standard  deviations  of  ZM  TEC  (left  plot)  calculated  from  the  CODE  data 
(magenta),  full  model  (green),  solar  part  of  the  model  (red)  and  temperature  part  of  the  model  (blue);  (right  plot)  the 
same  as  the  upper  plot  but  for  the  standard  deviations  of  DW  1  TEC. 
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In  order  to  compare  the  regression  model  with  the  TBC  data  variability  and  also  to  assess  the 
relative  effectiveness  of  the  temperature  and  F10.7  impact  on  the  TBC  variability  we  calculate  the 
STD  of  all  mentioned  components.  The  results  are  presented  in  Figure  36.  They  reveal  that:  (i)  the 
presented  linear  regression  model  describes  almost  half  of  the  real  variability  of  the  global  TBC  and 
(ii)  the  stratospheric  temperature  (i.e.  lower  atmospheric  forcing)  impact  is  almost  half  of  that  from 
the  F10.7  (i.e.  external  forcing  related  to  the  photo-ionization).  In  order  to  understand  better  the 
above  conclusions  we  remind  that  the  regression  model  includes  only  day-to-day  variability  of 
F10.7  and  temperature  (i.e.  wind  also).  It  is  worth  nothing  that  the  forcing  from  below  by  tides  and 
GWs  is  not  considered.  Having  in  mind  also  that  the  geomagnetic  activity  is  not  included  in  the 
model  then  it  is  understandable  why  the  model  can  describe  only  half  of  the  real  TBC  variability. 
The  contribution  of  the  temperature  dependence  to  the  TBC  variability  is  only  half  of  that  related  to 
F10.7  but  is  not  negligible  and  have  to  be  included  in  the  modelling  the  global  TBC. 


7.  Further  steps 

As  a  next  step  we  propose  to  build  a  new  global  TBC  model  for  short-term  prediction  where  the 
impact  of  stratospheric  temperature  is  included  as  well. 
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